Tag Info

Hot answers tagged

22

First off, your function is very simple without any hard number-crunching, so it will always be hard to get a large speedup for the compiled version. Secondly, your Parallelization option for Compile is useless because it doesn't do any parallelization this way. Let me give slightly changed versions of your examples and explain how you can achieve a large ...


21

Chunks of weak compositions Here is slightly modified version of algorithm used in Combinatorica`NextComposition converted to a LibraryFunction. Needs["CCompilerDriver`"] " #include \"WolframLibrary.h\" DLLEXPORT mint WolframLibrary_getVersion() { return WolframLibraryVersion; } DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) { ...


17

Here is a summary of comments (before @ciao's best answer above), with a change in notation. These functions calculate the number of partitions of n into exactly k distinct parts of size at most m. NumberOfWays000[n_, k_, m_] := Count[Map[Length,Map[DeleteDuplicates, IntegerPartitions[n,{k},Range[m]]]], k] NumberOfWays001[n_, k_, m_] := ...


16

This seems pretty quick, particularly on larger cases / larger k, e.g. 451, 29, 101 finishes in a few seconds on the loungebook. N.B. - I have not tested this exhaustively, just thrown together from ideas... If[Min[#3, #1 - Tr@Range@(#2 - 1)] < 0, 0, SeriesCoefficient[QPochhammer[-x y, x, Min[#3, #1 - Tr@Range@(#2 - 1)]], {x , ...


15

Here is a totally different approach based on the fact that successive products forming the generating function are due to multiplication by a binomial $1+t*z^j$. Form a matrix $v$ of zeros with $n+1$ rows and $k+1$ columns. Initialize the top left corner to 1. Iterate $v=v+w$ where $w$ is the matrix $v$ shifted down by $j$ rows and to the right by 1. The ...


10

I'm posting a whole new answer because I don't want to inherit any of the votes I received for my previous wrong answer. In formulating my new answer, I was aiming for correctness, simplicity, and reasonable (but not stellar) performance. Simplicity was achieved by taking a recursive approach, the clarity of which gives me confidence in the correctness of ...


9

The problem is that it reevaluates the sum every single time you call it, recomputing every 20^4 term again and again. You just need to compile the function CC2 so that it performs the summation only once. Using the code you have, it takes my machine about 6 seconds to compute a single data point: CC2[0.003] // AbsoluteTiming (* {6.069311, 1.49893} *) ...


7

This is faster for the given array: keepFirstK7[lst_List] := Delete[lst, Flatten[(GatherBy[Position[lst, #], Last] & /@ DeleteDuplicates[Flatten@lst])[[All, All, 2 ;;]], 2]] keepFirstK7@rArray $\ $ {{4, 2, 1, 0}, {2}, {0, 0, 3, 3, 3}, {}, {3, 4}} Another approach keepFirstK[lst_List] := Block[{sowFirst}, sowFirst[int_, count_] := ...


7

Another solution using memoization: keepFirst[lst_List] := Module[{isFirst}, isFirst[elem_, col_] := (isFirst[elem, col] = False; True); Pick[lst, MapIndexed[isFirst[#1, Last[#2]] &, lst, {2}]] ]


7

Well, to be honest, despite I've been using Mathematica for 3 years, I'm getting more and more confused about what's functional programming, but the following solution is at least more elegant and faster than yours: searchSpan2[knots_, u0_] := First@Ordering[UnitStep[u0 - knots], 1] - 1 NonzeroBasis2[p_, u_, u0_] := With[ {i = searchSpan2[u, u0], ...


6

To speed up the CAGDBezierSurface,I have two trials 1, Refactor the pure function AllBasis I don't know why the most efficient method to calculate the Benstein basis is the defintion of Benstein basis. AllBasisNew = Function[{deg, u0}, Bernstein[deg, #, u0] & /@ Range[0, deg]] AllBasisOld = Function[{deg, u0}, ...


6

Just a way. Sorts first and removes identical elements, then elements that are strict subsets of other elements... fun[lst_] := Module[{sort = DeleteDuplicates[Sort[Sort /@ lst]], sa, w}, sa = SparseArray@ Outer[Boole[SubsetQ[#1 /. w -> List, #2 /. w -> List]] &, w @@@ sort, w @@@ sort]; ReplacePart[sort, Thread[(Last /@ ...


6

You should use the SPICE toolkit (in C), which you can access from Mathematica with LibraryLink. The toolkit routines will handle all of the kernel access and coordinate transformations for you transparently and fast. In my opinion it would be both a waste of time and highly error prone to try to duplicate their functionality. Not to mention slower. The ...


5

A function for data fully-sorted: myDelete[data_] := Block[{revdata = Reverse@data, manip}, manip[l_] /; Length@l == 1 := l; manip[l_] := Block[{rest = Rest@l, subsetq, firstelem = First@l}, subsetq = SubsetQ[firstelem, #] & /@ rest; {firstelem, Sequence @@ manip[Pick[rest, subsetq, False]]}]; Reverse@manip[revdata]]; ...


5

NSolve with adequate precision works well $Version (* "10.2.0 for Mac OS X x86 (64-bit) (July 7, 2015)" *) eqns = Sin[z + Sin[z + Sin[z]]] == Cos[z + Cos[z + Cos[z]]] && -3 < Re[z] < 3 && -3 < Im[z] < 3; roots = NSolve[eqns, z, WorkingPrecision -> 20]; And @@ (eqns /. roots) (* True *) Note that there are a large ...


5

Compile f and use a memo-ized version of it Since it seems like NIntegrate decides to symbolically evaluate its argument first, I thought I'd force it not to by compiling the function f. This seems to make a significant difference: Clear[f, f1, g] g[x_] = Nest[f[x] + 1./# &, f[x], 500]; f1 = Compile[{x}, Sum[1/100 Erfc[-(x^2/k)], {k, 100}]]; ...


5

For a non-symmetric real matrix you can consider using LibraryLink to speed things up. It still won't be as fast as the Total/Tr answer, but it may be useful otherwise (call this C program SumUpperTriangle.c): #include "WolframLibrary.h" DLLEXPORT int SumUpperTriangle(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) { /* Variable ...


5

I found this question quite interesting, so I thought I would collect the answers contributed in comments for future reference and to have the question appear as answered in search. I generated a slightly bigger matrix to play with, and minimally modified the code to render it independent of the size of the matrix. I also compared timings of each method to ...


5

For the auxiliary fuction searchSpan[] , which came form the following algorithm of The NURBS Book. where $U=\{\underbrace {a,\cdots ,a}_{p+1},u_{p+1},\cdots u_{m-p-1},\underbrace {b,\cdots,b}_{p+1}\}, \quad n=m-p-1$ Here is a rule-based solution that I implemented according to The Toad's answer NonzeroBasis[{deg_, knots_}, u0_] := Module[{coeff, i, ...


5

You're after the Swinnerton-Dyer Polynomials. Take a look and compare with OEIS (which, BTW cites our friend Roman E. Maeder. Programming in Mathematica, Addison-Wesley, 1990, page 105): MinimalPolynomial[Sum[Sqrt[Prime[i]], {i, #}], x] & /@ Range@5 gives the same results shown in the OEIS page, but MinimalPolynomial can't calculate sixth term within ...


4

According to the solution of xzczd that dealing with the calculation of $$N_{i-p,p}(u_0),B_{i-p+1,p}(u_0), \cdots, N_{i,p}(u_0)$$ I mimic this strategy to calculate all the values of Berstern basis function of degree $n$ $$\color{blue}{B_{n,0}(u_0),B_{n,1}(u_0), \cdots, B{n,n}(u_0)}$$ AllBernsteinBasis[n_, u0_] := Nest[MovingAverage[ArrayPad[#, 1], {u0, ...


4

Evaluate the integral once and for all (cf. cdfc): cdfc[k_] = Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}]; TCJS[T_, k_] := A/T + c1[T]*d1*T + h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + (b1/T)*σ1* Sqrt[T + L1]*(PDF[NormalDistribution[0, 1], k] - k*cdfc[k]); EQ1[T_] := (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/ ...


4

At the end of this post you'll find the code for a small benchmark to compare LibraryLink with standard passing vs LibraryLink with MathLink based passing on two counts: function call overhead; use a function that does nothing, has no arguments, returns nothing (llNone and mlNone for LibraryLink and MathLink, respectively) passing real arrays; use a ...


3

Using: $posLenIdx = Association@MapIndexed[First@#2-> Length@#&, test]; $revItenIdx = Merge[Identity]@MapIndexed[Association@Thread[#-> First@#2]&, test]; and: testSet[set_List, {pos_}]:=Block[{biggerGroups, len = Length@set}, biggerGroups = Select[Tally@Flatten@Lookup[$revItenIdx, set, {}], ( #[[2(*qtd*)]] >= ...


3

I thought it interesting to ask where the roots determined by Bob Hanlon and Michael E2 lie in the complex plane. pts = Flatten[N[roots, 15] /. Rule[_, z_] -> ReIm[z], 1]; pts2 = Flatten[N[roots2, 15] /. Rule[_, z_] -> ReIm[z], 1]; As noted in their answers, the numbers of roots are 883 and 1251. One might suppose that the first list is a subset of ...


3

In V10, Solve works, too, and gives 1251 solutions. roots2 = Solve[eqns, z]; // AbsoluteTiming Length@roots2 Solve::incs: Warning: Solve was unable to prove that the solution set found is complete. >> (* {99.1951, Null} 1251 *) Maybe there are more, too. The timing is almost 6 times as long as BobHanlon's NSolve command on my computer. But ...


3

Here's a single-shot adaptation of my earlier (deleted since question changed) ideas. A bit over twice as fast as those for single-shot cases only, e.g. on c=5, 5000 sublists, 30 sublist length, p=0.001 done once. mutateList5[list_, c_, p_] := Module[{fl = Flatten@list, cr = Range@c, cl, lfl = Times @@ Dimensions@list, chg}, cl = Delete[cr, #] & /@ ...


3

Expanding @ciao's comment above: n = 3 10^6; n1 = 10000; v1 = Range[n]; o1 = Union@RandomInteger[{1, n}, n1]; o2 = Reverse[ o1/1000 // N]; (*Anything*) getV2[v1_, o1_, o2_] := Module[{f, v2}, v2 = 0 v1; v2[[o1]] = o1; (f[#1] = #2) & @@@ Transpose[{o1, o2}]; f[0] = o2[[1]]; v2 = f /@ FoldList[Max, v2] ] getV2[v1, o1, o2]; // Timing


3

My interpretation of this question is that the OP is asking for a zero-order interpolation of his data. There for I propose the following solution. First a proof of concept using about 50 points based on belisarius' data model. SeedRandom[42]; Module[{interp, n, n1, o1, o2}, n = 50; n1 = 15; o1 = {1, Sequence @@ Union@RandomInteger[{2, n - ...


3

I'd like to expand on the answer given by @ChipHurst in the comments. My hope was that one of these approximate tests would be much faster than PrimeQ; however, that was not the case. Perhaps someone can code these methods more efficiently than shown here. Fermat's Little Theorem: Prime p and any b such that GCD[b,p]=1 implies b^(p-1)=1, mod p. Use a prime ...



Only top voted, non community-wiki answers of a minimum length are eligible