Hot answers tagged

35

General comments First, if you plan to use multi-dimensional integrals it is better to test with multi-dimensional integrals not with one dimensional ones. One might think that the test in the question is an appropriate one if multi-dimensional integration is done by the integrator in a recursive manner. This seems to be case for scipy.integrate.nquad (see ...


33

I don't think it is possible for FullSimplify to assess how far it is from a meaningful reduction of complexity. At every stage parts of the expression may after some transformation cancel (or not). Sometimes complexity has to be increased before it can be lowered (as I showed in this answer). Simplification is the process of minimizing complexity, where ...


24

I can't take much credit for this answer--I hadn't even got version 10.2 installed until J. M. commented to me that these functions could be written efficiently in terms of the Hamming weight function. But, it is understandable that he doesn't want to write an answer using a smartphone. The definition of the built-in ThueMorse is: ThueMorse[n_Integer] := ...


24

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 ...


23

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) { ...


23

I've always wondered about the scalability of MathLink (now officially "Wolfram Symbolic Transfer Protocol"). This is the protocol used by Mathematica to communicate between the front end and the kernel, and the basis of the Parallel` package. It has quite low bandwidth and high latency relative to, for example, MPI libraries. I also wonder how many MathLink ...


22

I've decided to expand on my comment. Before I delve into the solution, let's all pause for a moment and marvel at the stereographic parametrization of a unit circle: $$\begin{pmatrix}\frac{1-t^2}{1+t^2}\\\frac{2t}{1+t^2}\end{pmatrix}$$ Sometimes also referred to as the Weierstrass substitution, it has often been used as a tool in the solution of algebraic ...


21

dat = {0.71, 0.685, 0.16, 0.82, 0.73, 0.44, 0.89, 0.02, 0.47, 0.65}; Module[{t = 0}, Split[dat, (t += #) <= 1 || (t = 0) &] ] {{0.71, 0.685}, {0.16, 0.82, 0.73}, {0.44, 0.89}, {0.02, 0.47, 0.65}} Credit to Simon Woods for getting me to think about using Or in applications like this. Performance I decided to make an attempt at a higher ...


20

We can take advantage of the fact that IntegerDigits is very fast when the base is large. But not too large: no bigger than $2^{63}-1$ on a 64-bit system or $2^{31}-1$ on a 32-bit one, because Mathematica's machine integers are signed. Additionally, non-power-of-two bases require more work to get the result than just partitioning a bit-string, and are ...


19

It definitely has something to do with the Interpolation function. Evaluating tempdata = Import["http://www.inrim.it/~magni/cm.dat.gz", "Table"]; cmfunc = Interpolation[tempdata] we get the warning Interpolation::udeg: Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will ...


19

SatisfiableQ has three methods: "BDD": converts the expression to a BDD (binary decision diagram), "SAT": uses the Minisat library, "TREE": a branch-and-bound method based on the expression tree. SatisfiabilityCount counts instances by converting the expression to a BDD, so its timing should be close to SatisfiableQ with the "BDD" method (counting ...


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_] := Total[Boole[...


17

I guess I should not have been surprised that there are actually many ways to estimate the Gaussian and mean curvature of a triangular mesh. I shall present here a slightly compacted implementation of the MDSB method, with a few of my own wrinkles added in. I should say outright that the current implementation has two weaknesses: the determination of the ...


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 , ...


16

You are right, it can be done in a fraction of second. One can explicitly construct an array of indexes blockArray[mat_] := SparseArray[ Tuples[Range@# - {1, 0, 0}].{Rest@#, {1, 0}, {0, 1}} &@Dimensions@mat -> Flatten@mat] Timings: matrices = RandomReal[1, {48, 128, 128}]; s1 = SparseArray@ ArrayFlatten@ReleaseHold@DiagonalMatrix[...


16

Vectorization will help a lot: a[x_?NumericQ] := N[Exp[-Abs[x]]]; x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]; A = Outer[a[#1 - #2] &, x, x]; // AbsoluteTiming (* {2.11988, Null} *) B = Exp[-Abs[x - #]] & /@ x; // AbsoluteTiming (* {0.016182, Null} *) A == B (* True *) Notice that I am doing arithmetic on vectors the size of x instead of scalars....


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 ...


15

The following is a much faster, but not optimal, recursive solution: pts = RandomReal[1, {10000, 2}]; f = Nearest[pts]; k[{}, r_] := r k[ptsaux_, r_: {}] := Module[{x = RandomChoice[ptsaux]}, k[Complement[ptsaux, f[x, {Infinity, .05}]], Append[r, x]]] ListPlot@k[pts] Some timings show this is two orders of magnitude faster ...


15

You are needlessly computing exactly the same BesselJ and BesselJZero function values over and over again. As an example, in the 100x100 case, BesselJZero[1 + l, 1] (l is a constant) is computed 551 times! You just need to compute each one once. The easy way to do that is to memoize: bj[n_?NumericQ, z_?NumericQ] := bj[n, z] = N[BesselJ[n, z]] bjz[n_?...


14

This is best I can do so far. The system is linear so LinearSolve is a natural thing to try. arrays = CoefficientArrays[eqs, vars] (* {SparseArray[< 2 >, {9}], SparseArray[< 35 >, {9, 9}]} *) solv1 = Thread[ vars -> LinearSolve[arrays[[2]], -arrays[[1]], Method -> "CofactorExpansion"]]; // AbsoluteTiming (* {0.28347, Null} ...


14

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 ...


14

pts = Partition[RandomReal[1, 10000], 2]; ListPlot[pts] Use SameTest option with Union pts2 = Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)]; Length[pts2] 326 ListPlot[pts2]


14

You can accomplish this using DeleteDuplicatesBy, by first taking your two input lists and making a matrix out of them, and then deleting the rows where the last element (the element that came from myNewList) is a duplicate. Then you transpose back and assign the sublists of the reduced matrix to the new lists you want. {myListDuplicatesDeleted, ...


14

I am sure you can easily install also Linux on it and then you could contact Vladyslav Shtabovenko, the current maintainer of FeynCalc (https://github.com/vsht) and ask him about hard problems in High Energy Physics he would like to benchmark on such a King-Kong machine. Either him or somebody else could also provide you with more complicated examples of ...


14

It took me a while, but the suggestion of @Michael E2 was quite helpful, and especially the post (Optimize inner loops). For those of you (like me) who are new to this style of programming in Mathematica there are a few things that helped in my particular example. In my slow version I was looping over all vertices in the mesh list. For example in ...


13

Outer is highly optimized for several built-in functions (Plus, Times, List). Therefore Exp@-Abs@Outer[Plus, #, -#] &@Range[-10, 10, 0.02]; // RepeatedTiming (* {0.025, Null} *) gives ~50x speedup over Outer[#1 - #2&, #, #] and ~15x speedup over Outer[Subtract, #, #]. Also is a bit faster then Kuba's Exp[-Abs[x - # & /@ x]].


13

Select[ IntegerPartitions[24, {8}, Range[5]], #.# == 86 & ] {{5, 5, 4, 2, 2, 2, 2, 2}, {5, 5, 3, 3, 3, 2, 2, 1}, {5, 4, 4, 4, 2, 2, 2, 1}, {5, 4, 4, 3, 3, 3, 1, 1}, {4, 4, 4, 4, 4, 2, 1, 1}} Slightly more general approach (in case where IntegerPartitions is not what we need): ClearAll[ar, a]; ar = Array[a, 8] ar /. Solve[Flatten@{ Tr[ar]...


13

That's not entirely true: a // Dimensions {2251, 101, 2} So it's a 2d array of paris, and each row is treated as separate set to plot, ListPlot is trying to style them differently, etc. That's what takes time probably. You can keep pairs and do a fast plot reducing it to simple list of pairs. With Flatten[#,1] or Catenate: ListPlot[ a ] // ...


13

Here's a million points processed in half a second: SeedRandom[0]; (* updated for reproducibility *) xyCoordinateCentreCircle = RandomReal[1, {1*^6, 3}]; Map[ Mean[#[[All, -1]]] &, BinLists[xyCoordinateCentreCircle, {0, 1, 1/10}, {0, 1, 1/10}, {0, 1, 1}], {3}] // AbsoluteTiming (* {0.513721, {{{0.506174}, {0.497757},..., {0.50284}}, {{...


12

I remember reading somewhere that everytime I use =, Mathematica copies an expression in the memory (which may be slow and inefficient). This is not quite true, as written here. Mathematica uses a copy-on-write behaviour, i.e. it will only create an actual copy of a datastructure if you modify it. Example: a = {1,2,3}; As this is evaluated, first the ...



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