Symbolic FAQ

Daniel Lichtblau
Wolfram Research, Inc.
danl@wolfram.com
Original version: June 1998
Most recent update: January 2005

Abstract

We will discuss several problems that appear from time to time as questions to tech-support or the Mathematica users' news group comp.soft-sys.math.mathematica. While not all are frequently asked, those that are not fall under the heading of "general interest." For example, we have some that involve essential computer-science ideas such as bit vector or binary tree or hash table construction and use.

I stole shamelessly from my own responses. For code not written by myself I give credit to  appropriate parties. All examples were run in a development version kernel.

\(Off[General::spell]; 
Off[General::spell1]; \)

Miscellaneous brief items

Q: How do I sort a list in such a way that I can obtain the position j occupied by the k-th element in the original ordering.

A: We augment the list with positions, sort the resulting ordered pairs, then Transpose to recover the sorted list of original elements and the corresponding positions from which they came.

sortWithPositions[ll_List] := {Sort[ll], Ordering[ll]}

ee = Table[Random[],  {10}];
sortWithPositions[ee]

Q: Why does Mathematica not remove the common factor in

Cancel[(3. *a^3  +  a)/(2*a^2  +  a)]

A: Algebraic operations are very difficult to implement when working with inexact numbers. While the above is a glaring case of a problem that should give an obvious simplification (and may in a future release of Mathematica), this is generally not well understood. If you want algebraic functionality you generally need to use exact coefficients.

Q: How can I find the rank of a matrix?

A: (If you asked me directly, the actual answer might be: "Depends on how you lost it"). The rank is the number of pivots in the (row) reduced echelon form.

\(rank[mat_?MatrixQ]\  := \ 
	
      Length[Select[RowReduce[mat], \ Max[#] =!= 0&]]\)

rank[{{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}]

If you have a numeric matrix there may be considerations of numeric stability. In this case one might use SingularValues, if the matrix is square (or augment it to make it square if not), or perhaps use RowReduce with a zero test that uses, say, Chop.

Q: I have a 50 digit number and I want to pull out any factors that are reasonably quick to find. Does Mathematica provide a way to do this?

Yes. Simply use FactorInteger[...,  FactorComplete->False]. In the example below we will generate a number with some prime factors of 12-13 digits, then see how quickly one may pull out the smaller digits.

aa = 
      Join[Table[j, {j, (10^12) + 39, (10^12) + 43, 2}], 
        Table[j, {j, (10^13) + 51, (10^13) + 57, 2}]]PrimeQ[aa]bb = Times@@aa;

FactorInteger[bb,  FactorComplete -> False]  //  Timing

Note that there may be more easy factors to pull out from the ones already obtained.

FactorInteger[
        7607197928576694595777548789277900635696847387955303994636410490017,  FactorComplete -> False]  //  Timing

FactorInteger[257265294119113471,   FactorComplete -> False]  //  
      Timing

FactorInteger[29569468181178657565420783168444872229695315935327,  
        FactorComplete -> False]  //  Timing

At this point one can check that all factors extracted thus far, except the last one above (which has 42 digits), are prime. It takes about half a minute to factor that one on my pentium pro.

Q: I have a list of starting values that I want to use in FindRoot. How might I do this without retyping them one by one, given that FindRoot does not want a list of these but rather separate entries for each?

A: You can apply Sequence to the list in order to splice in the elements as separate arguments to some function. In the example below, the Evaluate wrapper is needed because FindRoot has the attribute HoldAll.

\(startlist\  = \ {{x1, \ 1}, \ {x2, \ 1}, \ {x3, \ 1}}; 

    eqns\  = \ {x1^2 - 2*x1*x2 - x2*x3 + 3 == 0, \ x1*x3 == x2, \ 
        x1 + x2 - x3^2 == 3}; 

    FindRoot[eqns, \ Evaluate[Apply[Sequence, \ startlist]]]\)

Q: How do I decide when to use Solve, when Reduce, when LinearSolve, when FindRoot, when ....?

A: This and other Solve-related issues are addressed in the notebook EquationSolving.nb

Q: How do I "collect" with respect to arbitrary powers of a variable?

A: One way to perform the Collect is to allow a pattern variable for the exponent. This will in effect create a separate "variable" for each distinct power, and this may suffice to do what you want. For  example:

\(expr\  = \ 
      Expand[Sum[\((a*n\  + b*n\  - c)\)*x^\((n + 0.12/n)\), {n, 1, 4}]]; 

    Collect[expr, \ x^_]\)

Q: How do I "collect" with respect to real and imaginary parts separately?

A: Here is a simple method.

complexCollect[expr_, var_] := 
      Module[{myI, e2, 
          v2}, {e2, v2} = {expr, var} /. 
            Complex[a_, b_] -> a + myI b; 
        Collect[e2, v2] /. myI -> I]

complexCollect[(a + I) b  +  c I ,  I]

Q: How might I obtain the Cartesian form of an expression?

A: The way to do this at present is via ComplexExpand. Unless told otherwise it assumes all variables are real-valued.

\(expr\  = \ \((1 + 2  I)\)*Exp[I*x]; 

    ComplexExpand[{Re[expr], \ Im[expr]}]\)

Here is another example, one that apparently does nothing useful.

ComplexExpand[Abs[x + I*y]]

We opted to allow Abs, etc. be used in  results except when explicitly excluded by not appearing in a given TargetFunctions list. So returning Abs[x+I*y] is okay since the job of ComplexExpand is to separate into real and imaginary parts and Abs[...] is indeed explicitly real. To get a different result one can do as below.

ComplexExpand[Abs[x + I*y],  TargetFunctions -> {Re, Im}]

Examples of polynomial manipulation

Q: I want to expand multinomials directly, that is, without recourse to Mathematica's (excellent, mind you) Expand function. What are some ways in which this might be done?

Here are a few ways to do expansion to integer powers, coded in Mathematica. The first is straightforward multiplication/expansion in a loop up to the degree desired. The second is the divide-degree-by-two approach, working recursively. The third is similar, but works from the bottom up, building what is called a the power-tree. The fourth is brutal multinomial expansion. Another way to do such, not shown, might be via integer partitions.

I include two examples. Of course these do not begin to cover the vast realm of "reasonable" inputs one might want to test. Note that several of these routines use the built-in Mathematica Expand function. In all cases I do so either to expand a square, or to expand the product of two (typically) different polynomials. One could replace these by more elementary code that avoids such higher-level functionality, but I doubt the results would change appreciably in regards to performance.

simpleExpand[deg_Integer, poly_] := 
      Module[{p = poly, result = poly}, 
        Do[result = Expand[result poly], {deg - 1}]; result]

recursiveExpand[deg_Integer, poly_] := 
      Module[{result}, If[deg === 1, Return[poly]]; 
        If[deg === 2, Return[Expand[(poly^2)]]]; 
        result = recursiveExpand[Floor[(deg/2)], poly]; 
        result = Expand[(result^2)]; 
        If[OddQ[deg], result = Expand[result poly]]; result]

powerTreeExpand[deg_Integer, poly_] := 
      Module[{psquare = poly, j = deg, result = 1}, 
        While[Quotient[j, 2] =!= 0, 
          If[Mod[j, 2] === 1, result = Expand[result psquare]]; 
          psquare = Expand[(psquare^2)]; j = Quotient[j, 2]; ]; 
        Expand[result psquare]]

multinomialExpand[n_Integer, x_List] := 
      Module[{len = Length[x], jj, indices, expons, iterators}, 
        indices = Array[jj, {len - 1}]; 
        expons = {Sequence@@indices, n - Plus@@indices}; 
        iterators = 
          Table[{indices[[j]], 0, 
              n - Plus@@Take[indices, j - 1]}, {j, len - 1}]; 
        Plus@@Flatten[
            Table[Multinomial[Sequence@@expons] Times@@(x^expons), 
              Evaluate[Sequence@@iterators]]]]

Timing[smallex1 = simpleExpand[10, k + r + x + k r]; ]Timing[smallex2 = recursiveExpand[10, k + r + x + k r]; ]Timing[smallex3 = powerTreeExpand[10, k + r + x + k r]; ]Timing[smallex4 = multinomialExpand[10, {k, r, x, k r}]; ]smallex1 === smallex2 === smallex3 === smallex4 === 
      Expand[(((k + r + x + k r))^10)]

Timing[bigex1 = simpleExpand[20, k + r + x + k r]; ]Timing[bigex2 = recursiveExpand[20, k + r + x + k r]; ]Timing[bigex3 = powerTreeExpand[20, k + r + x + k r]; ]Timing[bigex4 = multinomialExpand[20, {k, r, x, k r}]; ]bigex1 === bigex2 === bigex3 === bigex4 === 
      Expand[(((k + r + x + k r))^20)]

Q: I have a very large polynomial in Mathematica that has several "small" variables. How might I discard all but lowest order terms, making no assumptions regarding the relative sizes of the small variables? For example given 40 (x^3) y z + 11 x (y^2) z - 18 x (y^3) z + 
          28 (x^3) y (z^2) with all variables considered "small," I would like to get as a result just the first two terms 40 (x^3) y z + 11 x (y^2) z because the third term is dominated by the second and the fourth term is dominated by the first.

A: You want to keep only those terms whose power-products are not divisible by those of other terms. These can by found by getting a Gröbner basis of the set of monomials.

polyxyz = 
      40 (x^3) y z + 11 x (y^2) z - 18 x (y^3) z + 28 (x^3) y (z^2); 
    ff = GroebnerBasis[List @@ polyxyz]
gg = Coefficient[polyxyz, ff]
gg . ff

Q: How do I "collect" an expression algebraically with respect to products or powers? For example, I want to collect terms in a (x^2) + a b x y with respect to a*x so that I get a*x*(x+b*y).

A: We can then do the job with PolynomialReduce. We first use it to replace, algebraically, all occurences of our products and powers with new variables. We then collect with respect to the new variables (using Collect).  Last we again use PolynomialReduce, now to rewrite the new variables in terms of the original ones. When we collect we will encapsulate the coefficients so that the latter PolynomialReduce does not expand them. This cannot be done by setting the CoefficientDomain option to RationalFunctions because the coefficients may involve the same variables as are used in the products and powers of interest; instead we must use a dummy head.

algebraicCollect[expr_, terms_List] := 
      Module[{temp, newexpr, newvars, repvars, nv, redpolys, vars, result, 
          plus}, newexpr = Prepend[terms, expr]; 
        newvars = Variables[newexpr]; repvars = Array[nv, Length[terms]]; 
        vars = Join[newvars, repvars]; 
redpolys = 
          terms - repvars; 
        result = Last[PolynomialReduce[expr, redpolys, vars]]; 
        result = Collect[result, repvars, plus]; 
        Last[PolynomialReduce[result, redpolys, 
              Reverse[vars]]] /. plus -> 
            Plus]
algebraicCollect[expr_, term_] := 
      algebraicCollect[expr, {term}]

Here are some simple examples.

algebraicCollect[a (x^2) + a b x y + c (y^2) + a y + b x, a x]
      algebraicCollect[
      (x^5) + a (x^3) + a b (x^4) y + c (x^2) (y^2) + a x y + b x, (x^2)]

Q: I have a polynomial in the variables x,y,z. I would like to replace algebraically all occurrences of (y^2) by z - (x^2) - 2. If I call PolynomialReduce directly it does not always work? What is the correct way to do this?

A :  Here is a simple example .

poly = 
      -(x^3) + 9 x (y^2) + 2 (y^3) + x (y^3) + 6 (x^2) z + (x^3) z - 
        3 y (z^2) + 3 (z^3) + y (z^3); Last[PolynomialReduce[poly, (y^2) - (z - (x^2) - 2)]]

Apparently this did not work. Why? Because we failed to specify a term order in which (y^2)would be "larger" than the rest of the reducing polynomial. To do this it suffices simply to list y as the variable of interest.

Last[PolynomialReduce[poly, (y^2) - (z - (x^2) - 2), y]]

The general functionality responsible for algebraic "replacement" is GroebnerBasis/PolynomialReduce. One sets up canonical replacement polynomials with GroebnerBasis, then uses these as the second argument to PolynomialReduce (making certain to use the same term ordering). For example, suppose you want to rewrite poly above in such a way as to minimize total degrees, given the "rules" {(x^3) + x (y^2) + (z^2) - 2, x z + 2 (y^2) - 3 x y - 1, 
        (x^2) + 2 y z - 3}.

gb = 
        GroebnerBasis[
          {(x^3) + x (y^2) + (z^2) - 2, x z + 2 (y^2) - 3 x y - 1, 
            (x^2) + 2 y z - 3}, {x, y, z}, 
          MonomialOrder -> DegreeLexicographic];

Last[
      PolynomialReduce[poly,  gb,  {x, y, z},  
        MonomialOrder -> DegreeLexicographic]]

Q: Suppose I want to do my algebraic replacement in such a way that lower degrees get replaced by higher. How might I do this?

A: So now you want to rewrite in such a way as to maximize total degree. This can be done via an interesting tactic. First "homogenize" both input and reducing polynomials with respect to a new variable. Compute a Gröbner basis with a total-degree term order, such that the new variable is lexically "largest". When we later dehomogenize this will have the effect of making lower-degree monomials reduce in terms of higher-degree monomials. The code to homogenize shown below assumes that the input polynomial is expanded with respect to the variables of interest.

homogenize[poly_, vars_, 
          new_] /; (Head[poly] =!= Plus  &&  Head[poly] =!= List) := 
      poly
homogenize[poly_, vars_, new_] := 
      Module[{degfunc, totdeg, j}, 
        degfunc = 
          Apply[Plus, 
            Map[Function[x, Exponent[#, x] &], 
              vars]]; 
degfunc = 
          Distribute[degfunc, Function, Plus]; 
totdeg = 
          Max[Map[degfunc, Apply[List, poly]]]; 
Apply[
          Plus, Table[
            poly[[j]]*new^(totdeg - degfunc[poly[[j]]]), {j, 
              Length[poly]}]]]
homogenize[polys_List, vars_,  v_] := 
      Map[homogenize[#, vars,  v] &, polys]

redset = 
      homogenize[{(x^3) + x (y^2) + (z^2) - 2, x z + 2 (y^2) - 3 x y - 1, 
          (x^2) + 2 y z - 3}, {x, y, z}, w]

hgb  =  
      GroebnerBasis[redset,  {w, x, y, z},  
        MonomialOrder -> DegreeReverseLexicographic];

Last[
        PolynomialReduce[homogenize[poly, {x, y, z}, w],  hgb, 
           {w, x, y, z},  MonomialOrder -> DegreeLexicographic]]  /. 
       w -> 1

Some noncommutative algebraic manipulation

Q: How might one form and use the operator p(t,d/dt):=(t*d/dt-a)°(t*d/dt-b)°(t*d/dt-c) ?

A: Note that multiplication of t and d/dt do not commute. It is strongly recommended that one NOT change the attributes on Times to handle this. It is perhaps best just to define a differential operator function. This is the approach we take below. We will use the built-in NonCommutativeMultiply, although we could just as well make up our own operator. In fact we do exactly this to handle Lie bracket via bracket[]. One could instead enforce different commutator rules.

Clear[differentialOperate]
differentialOperate[a_, expr_] /; FreeQ[a, D] := 
      a*expr
differentialOperate[L1_ + L2_, expr_] := 
      differentialOperate[L1, expr] + 
        differentialOperate[L2, expr]
differentialOperate[a_*L_, expr_] /; FreeQ[a, D] := 
      a*differentialOperate[L, expr]
differentialOperate[a : HoldPattern[D[__] &], expr_] := 
      a[expr]
differentialOperate[L1__ ** L2_, expr_] := 
      Expand[differentialOperate[L1, 
          differentialOperate[L2, expr]]]
commutator[L1_, L2_] := L1 ** L2 - L2 ** L1
differentialOperate[L1_^n_Integer, expr_] /; n > 1 := 
      Nest[Expand[differentialOperate[L1, #]] &, expr, n]

We now may form a differential operator and apply it to an expression. Returning to the  example specified in the question, we obtain the following.

ddt  =  (D[#, t] &);
diffop = (t ** ddt - a) ** (t ** ddt - b) ** (t ** ddt - 
            a); differentialOperate[diffop, f[t]]

L0 = x;
ddx  =  (D[#, x] &);
L1 = a*x^(3/2)*ddx + b*x^2*ddx^2;
differentialOperate[commutator[L0, L1], f[t, x]]

e1 = differentialOperate[L1, f[t, x]]
e2 = differentialOperate[L1, e1]
Expand[differentialOperate[L1^2, f[t, x]] - e2]

L0 = Function[{x}, x];
L1 = 
        Function[{x}, 
          a*x^(3/2)*(D[#, x] &) + 
            b*x^2*(D[#, x] &)^2];

differentialOperate[commutator[L0[t], L1[t]], f[s, t]]

e1 = differentialOperate[L1[t], f[t, x]]
e2 = differentialOperate[L1[t], e1]
Expand[differentialOperate[L1[t]^2, f[t, x]] - e2]

Q: (from an old benchmark test suite): How might one simplify a product or power in a noncommutative algebra that has commutator rules? Specifically, we have variables of the
form s[j], j an integer in the range from 0 to 9, and we have the relations s[i]^2 == 1
 and s[i]*s[j]==-s[j]*s[i] for i not equal to j. We want to expand and simplify the fifth power of (s[0]+...+s[9]).

A: A quick and dirty method is shown below. We make times itself noncommutative and as we expand we employ the canonicalizing rules to simplify the intermediate results.

ClearAttributes[Times, Orderless]rules = {s[i_Integer]^2 :> 1, 
        a___ s[i_Integer] s[j_Integer] b___ /; j > i :> 
          -a s[j] s[i] b}expr = (∑_(j = 0)^9) s[j]Fold[Expand[#1 #2] //. rules&, expr, Table[expr, {4}]]

Again, this method is NOT recommended because clearing attributes in basic functions is not safe. One rarely knows in advance whether internal code will then go wrong because it relies on basic functions such as Times. And what happens to subsequent comptations if, once done, you forget to do SetAttributes[Times,Orderless]?

A better way is to define our own concommutative multiplication operator. We will give it several basic attributes and also tell it how to canonicalize according to the relations above. The attributes we want it to have are
(i) The empty product is 1.
(ii) The product of a singleton is that singleton.
(iii) The product is associative.
(iv) We will automatically distribute over Plus.
(v) Integer factors are to be pulled out and multiplied using ordinary
multiplication.
(vi) We will have it recognize the relation that s[i]^2 is one.
(vii) We will have it canonicalize the ordered product of s[i] and s[j],
if j is greater than i, to become the negative of the ordered product of
s[j] and s[i].
We implement these quite simply below.

ncTimes[]  :=  1ncTimes[a_]  :=  ancTimes[a___, ncTimes[b_, c__], d___]  :=  ncTimes[a, b, c, d]ncTimes[a___, x_ + y_, b___]  :=  
      ncTimes[a, x, b]  +  ncTimes[a, y, b]ncTimes[a___, i_Integer*c_, b___]  :=  i*ncTimes[a, c, b]ncTimes[a___, i_Integer, b___]  :=  i*ncTimes[a, b]ncTimes[a___, s[i_Integer], s[i_], b___]  :=  ncTimes[a, b]\(ncTimes[a___, s[i_Integer], s[j_Integer], b___]\  /; \ j > i\  := 
\ \ 
      \(-ncTimes[a, s[j], s[i], b]\)\)

Now let us redo the example (all the while feeling superior because we did not muck with attributes of basic arithmetic functions).

expr  =  Sum[s[j],  {j,  0,  9}]Fold[ncTimes[#1,  #2]&,  expr,  Table[expr,  {4}]]

Q: How might one implement algebraic manipulations involving scalars and anticommuting variables? We will designate all variables from the latter category with a special head e.g. "avar".

A: We will show code below that repairs some deficiencies in earlier versions that appeared in MathGroup and The Mathematica Journal.

myTimes[a___, avar[x_], b___, avar[x_], c___] := 0
myTimes[avar[x_], aa : avar[_] .. ] /; 
        Not[OrderedQ[{x, aa} /. (avar[f_]) :> 
                f]] := (Signature[{x, aa} /. (avar[f_]) :> 
                f]*myTimes[
            Sequence @@ 
              Map[avar, 
                Sort[{x, aa} /. (avar[f_]) :> 
                      f]]])
myTimes[a_] /; Head[a] =!= Times := a
myTimes[x___, y_Plus, z___] := myTimes[x, #, z] & /@ y
myTimes[x___, y_Times, z___] := 
      myTimes[x, Apply[Sequence, y], z]
myTimes[x___, y_myTimes, z___] := 
      myTimes[x, Apply[Sequence, y], z]
myTimes[] = 1;
myTimes[a_, b_] /; Head[a] =!= avar || Head[b] =!= avar := 
      a*b
myTimes[vars__] /; Length[{vars}] > 2 := 
      Module[{vl = {vars}, antis, coms}, 
        antis = Cases[vl, avar[_]]; 
coms = 
          DeleteCases[vl, _avar]; 
Apply[Times, coms]*
            Apply[myTimes, antis] /; Length[coms] > 0]

myTimes[avar[x], a + b + c + d, z, avar[u]]

ncDet[{{a_}}] := a
ncDet[mat_?MatrixQ] /; 
        Length[mat] === Length[mat[[1]]] := ncDet[mat] = 
        With[{len = Length[mat], mat1 = Drop[mat, {1}]}, 
          Sum[(-1)^(j + 1)*
              myTimes[mat[[1, j]], ncDet[Drop[mat1, None, {j}]]], {j, 
              len}]]

m1 = {{avar[1], x, avar[2]}, {x2, avar[1] + 2*avar[2], 
            3*avar[1] - avar[2]}, {avar[1] - avar[2], 
            x3, -3*avar[1] + avar[2]}};
ncDet[m1]

Shuffling a large set efficiently

Q: How might I efficiently shuffle a large set in such a way that all permutations are equally likely?

A: I gave some mediocre answers to this in comp.soft-sys.math.mathematica a few months back. Let me try to improve on this. Say the set has length n. A simple method is to iterate over the elements from first to nth.. For element j, choose a random index k between j and  n inclusive. Swap the jth element with the kth.
Exercise 1: Why does this produce a shuffle such that all permutations are equally likely?
Exercise 2: Suppose that for element j we choose random k between 1 and n inclusive. Why does that NOT give equally probable shuffles? (Hint: consider how many outcomes there would be for n equal to 3 and how many permutations there are of 3 elements).

\(shuffle1[n_Integer]\  := \ 
      Module[{res = Range[n], \ tmp}, 
	
        Do[
		rand\  = \ Random[Integer, \ {j, n}]; 
		
          tmp\  = \ res[\([j]\)]; 
		res[\([j]\)]\  = \ res[\([rand]\)]; 

          		res[\([rand]\)]\  = \ tmp, 
		{j, 1, n}]; 
	res
	]\)

Timing[deck10ka  =  shuffle1[10^4]; ]Timing[deck20ka  =  shuffle1[2*10^4]; ]

In the currently available versions of Mathematica this is considerably slower. In fact one can check that it has quadratic rather than linear complexity. This is because the Part evaluations must reevaluate the entire list (or at least check that this is not necessary). Why is this not true of our development version? Because we work with a packed tensor of machine numbers (of type integer). So long as we set parts of the list to be integers this does not change, the entire list is a raw data type for which reevaluation is instantaneous.

How might one remove the need for reevaluation in the current version? Michael Trott suggested the following adroit use of Hold.

\(shuffle2[n_Integer]\  := \ 

      Module[{res = {Range[n]}, \ tmp}, 
\ \ \ \ \ \ \ \ 
        res\  = \ Apply[Hold, res]; 
\ \ \ \ \ \ \ \ 
        Do[rand\  = \ Random[Integer, \ {j, n}]; 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ tmp\  = \ res[\([1, j]\)]; 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          res[\([1, j]\)]\  = \ res[\([1, rand]\)]; 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ res[\([1, rand]\)]\  = \ tmp, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {j, 1, n}]; 
\ \ \ \ \ \ \ \ 
        res[\([1]\)]
\ \ \ \ \ \ \ \ ]\)

Timing[deck10kb  =  shuffle2[10^4]; ]Timing[deck20kb  =  shuffle2[2*10^4]; ]

One might notice that even though we have essentially linear complexity now this still seems to be slow, around two seconds to handle 10000 elements. If you comment out various parts of the module to find the culprit, you learn that all statements, Random, Part, and SetPart, are implicated. In other words, we are simply putting the core evaluator loop through alot of work. We can improve this immensely using Compile.

\(shuffle3\  = \ 
        Compile[{{n, _Integer}}, 
	
          Module[{res = Range[n], \ tmp, \ rand}, 
	
            Do[
		rand\  = \ Random[Integer, \ {j, n}]; 
		
              tmp\  = \ res[\([j]\)]; 
		
              res[\([j]\)]\  = \ res[\([rand]\)]; 
		
              res[\([rand]\)]\  = \ tmp, 
		{j, 1, n}]; 
	res
	]]; 

      shuffle3 // InputForm\) (* make sure it all got "compiled" *)

shuffle3[10^4];  // Timing

Actually we can improve this still further. When we obtain small random integers (fewer than 31 bits) we use a cellular automaton algorithm that generates one bit at a time. We might instead generate a random machine double and normalize so that it lies in the appropriate range, then take the integer floor.

\(\(shuffle4\  = \ 
      Compile[{{n, _Integer}}, 
	
        Module[{res = Range[n], \ tmp, \ rand}, 
	
          Do[
		rand\  = \ j\  + \ Floor[Random[]*\((n + 1 - j)\)]; 
		
            tmp\  = \ res[\([j]\)]; 
		res[\([j]\)]\  = \ res[\([rand]\)]; 
            
		res[\([rand]\)]\  = \ tmp, 
		{j, 1, n}]; 
	res
	]]; 
    \)\)shuffle4  //  InputForm

shuffle4[10^4];  // Timing

We'll verify that the complexity is O(n).

shuffle4[10^5];   //  Timing

shuffle4[10^6];   //  Timing

Of course if one only wants to pick a subset of k elements at random (say, dealing 20000 cards from a deck of 100000), the complexity is O(k); we do not need to pre-shuffle the entire deck at once. Minor alterations to the above routines will allow us to deal a subset of elements, then deal more, etc. Just like poker.

Bitwise operations and dense representation and manipulation of sets

Q: How might I obtain bitwise logical operations on nonnegative integers?

A: This functionality is built into our development version. For the present you might instead do as below.

bitop[opcode[aox_], a_, b_, c__] := 
      bitop[opcode[aox], bitop[opcode[aox], a, b], c]\(bitop[opcode[aox_], a_Integer, b_Integer] /; 
        \((a == 0 || a == 1)\) && \((b == 0 || b == 1)\) := 
      Switch[aox, 
and, a\ b, 
or, a + b - a\ b, 
xor, Mod[a + b, 2]]\)bitop[opcode[aox_], a_List, b_List] := 
      Module[{a1, b1, lena = Length[a], lenb = Length[b], len}, 
        len = Max[lena, lenb]; 
        a1 = If[lena < len, Join[Table[0, {len - lena}], a], a]; 
        b1 = If[lenb < len, Join[Table[0, {len - lenb}], b], b]; 
        Inner[bitop[opcode[aox], #1, #2]&, a1, b1, List]]bitop[opcode[aox_], a_Integer, b_Integer] := 
      FromDigits[
        bitop[opcode[aox], IntegerDigits[a, 2], IntegerDigits[b, 2]], 2]

Here are several examples.

{a1,  a2}  =  Table[Random[Integer,  10^10],  {2}]

bitop[opcode[and],  a1, a2]BitAnd[a1, a2]

bitop[opcode[or],  a1, a2]BitOr[a1, a2]

bitop[opcode[xor],  a1,  a2]BitXor[a1, a2]

Q: How might I represent elements of the power set of a set of n elements? I will need to perform logical operations such as Union and Intersection.

For such purposes one may represent elemets of the power set by integers in the range from 0 to (2^n) - 1. A value m in this range can represent a given subset by its binary string, where 1's stand for those elements in the subset.

For example, say we have a set of 100 elements. We will use integers in the range from 0 to (2^100) - 1 to represent subsets, and take a random subset of, say, 128 elements in the power set (that is, subset whose elements are subsets of the original set). Actually, there is nonzero probability that my subset contains stricly fewer than 128 elements, but this is still okay for the examples below.

k = 100; len = 128; 
    subs = Union[Table[Random[Integer, {0, (2^k) - 1}], {len}]];

A reasonable question might be: "How many elements of subs contain the third element of the original set?" To answer this we represent the elements of the original set as "bit masks", that is, integer powers of two.

masks = Table[(2^m), {m, 0, k - 1}]; k =.

So the third element of the set is masks[[3]], which is 2^2, or 4. We form the function that takes an integer and computes its bit-wise AND with masks[[3]]. We Map this function over the subset list subs, then count the number of elements in the result that are positive. The bit-wise operation has masked out all but the third binary digit of each member of subs, so the positive values resulting (all will be 4) are from those member that have a third binary digit set.

\(func\  = \ BitAnd[masks[\([3]\)], #]&; 


    Count[Map[func, \ subs], \ _?Positive]\)

Okay, easy enough. How about "How many members of subs contain the sixth and eleventh elements of the original set, but not the fourteenth?" (one might expect on average that this would be about an eighth of the number of elements in subs). Same idea.

func2  =  
        BitAnd[masks[[6]] + masks[[11]] + 
              masks[[14]],  #] &;
\(\(crit\  = \ \((# == 
              masks[\([6]\)] + masks[\([11]\)])\) &;\)
\)
Count[Map[func2,  subs],  _?crit]

What about "How many elements contain at least two from among the fourth, tenth, twenty-eighth, and fifty-first elements of the original set?" As before, we first mask appropriately.

func3  =  
      BitAnd[masks[[4]] + masks[[10]] + masks[[28]] + 
            masks[[51]], #]&;

Now we simply Count all those masked values whose DigitCount is at least two (by default the new function DigitCount will look at binary digits).

Count[Map[func3,  subs],  _?(DigitCount[#] >= 2&)]

Linear Algebra in a Galois field

Q: How can I find the row echelon form of

mat = {{a 2 x + 1, 2 a (x^2) + (a + 6) x + 5, 
          ((a^2) + a + 6) x + (3 a + 4)}, {
          (4 a + 2) (x^2) + (a + 2) x + 3, 
          ((a^2) + 3 a + 3) (x^2) + (a + 4) x + a, x + 1}};

in a field of characteristic 7 extended by an algebraic 
element a that satisfies

defpoly = (a^3) + 5 (a^2) + 3 a + 4;

A: Built-in functions will only do row reduction modulo primes for integer 
matrices, hence cannot support algebraic or transcendental 
elements (a and x respectively in the matrix above). We can 
do the job with GroebnerBasis instead.

newvars = {z1, z2, z3};gb = 
      GroebnerBasis[Append[mat . newvars, defpoly], Append[newvars, a], 
        CoefficientDomain -> RationalFunctions, Modulus -> 7]

The first element in the 
Gröbner basis will be the defining polynomial, so we drop it.

To return to matrix form we use Outer[D,polys,vars]. It is then a 
simple matter to get the rows ordered correctly. Sort will put 
rows beginning with more zeroes in front of those with less, so 
just reverse this.

redmat = Reverse[Sort[Outer[D, Drop[gb, 1], newvars]]]

Finally one might want to normalize rows so that pivots are 1.

normalize[vec_?VectorQ] := 
      With[{piv = Select[vec, # =!= 0 &, 1]}, 
        If[Length[piv] == 1, vec/piv[[1]], vec]]Map[normalize, redmat]

Minors vs. cofactors

Q: I want to get the cofactors of a matrix. The Minors command gives me what I want, but in a different ordering and with signs not alternating as needed. How can I get cofactors from Minors? (This is worded a bit euphemistically, if you will. We generally get this question as a bug report. But it is not a bug, as will be explained in the answer section).

A: First, note that minors are returned in lexicographic ordering with respect to rows and columns used to form them. In other words (as stated in the documentation), each element in the result corresponds to taking rows and columns with particular lists of positions. The ordering of the elements is such that reading across or down the final matrix the successive lists of positions appear in lexicographic order.

Mathematica will give minors for any size k less than or equal to the lesser of row and column dimension. Only in the special case where row dimensionis equal to column dimension and k is one less than this dimension does there seem to be confusion. In this case one obtains the cofactors rather simply from the minors by alternating signs and by reversing order of both rows and columns. This corresponds to a reverse-lexicographic ordering based on rows and columns used to form the cofactor (alternatively, one can describe it as lexicographic ordering based on the rows and columns not used).

In general we can relate cofactors to minors as follows. Take the  minor formed from the submatrix comprised by rows {r1,...,rk} and columns {c1,...,ck}. Then we can also form a (n-k)x(n-k) submatrix from the remaining rows and columns. Call the determinant of this submatrix d. Then the cofactor corresponding to rows {r1,...,rk} and columns {c1,...,ck} is (-1)^(r1+...+rk+c1+...+ck)*d. A well known special case is when k is one. Then the cofactors are formed from (n-1)x(n-1) submatrices, and the cofactor corresponding to the entry in row i and column j is the determinant of the matrix formed by removing row i and column j from the original matrix. Identities involving minors and cofactors are given by the Laplace expansion theorem. Perhaps best known of these states that, for k equal to one, the sum of products of 1x1 minors (that is, individual elements) in a given row times their corresponding cofactors gives the determinant of the matrix, while the sum of products of such minors by the cofactors from a different row gives zero. This is exactly the result used in elementary linear algebra to use the determinant and dimension n-1 cofactors to invert a matrix. The Laplace expansion result perhaps explains the origin of the term "cofactor."

Here is code to transform minors to cofactors.

rowSignMatrix[len_] := 
      Table[(-1)^(j), {j, len}, {k, len}]colSignMatrix[len_] := 
      Table[(-1)^(k), {j, len}, {k, len}]kSignMatrix[len_, k_] := 
      Module[{rs = rowSignMatrix[len], cs = colSignMatrix[len], ks1, ks2}, 
        ks1 = Minors[rs, k, Apply[Times, Transpose[#][[1]]]&]; 
        
ks2 = Minors[cs, k, Apply[Times, #[[1]]]&]; 
        
ks1*ks2]cofactors[mat_?MatrixQ, k_] := 
      Module[{len = Length[mat], ks, mn}, 
        (ks = kSignMatrix[len, len - k]; 
            
mn = Minors[mat, len - k]; 
            
mn = ks*mn; 
            
Reverse[Map[Reverse, mn]]) /; 
          len ≥ k]cofactors[mat_]  :=  cofactors[mat,  Length[mat] - 1]

Here is a simple example.

SeedRandom[1111]mat = Table[Random[Integer, 5], {3}, {3}];mn  =  Minors[mat]cf  =  cofactors[mat](mn . Transpose[cf])/Det[mat]

The dot product of minors by cofactors matrices must give the matrix determinant times an identity matrix of the appropriate dimension (this is the general statement of the Laplace expansion). We check this for a random 6x6 example.

SeedRandom[1111]mat = Table[Random[Integer, 5], {6}, {6}];checkLaplaceExpansion[mat_] := 
      Module[{m, c, det = Det[mat]}, 
        Table[m = Minors[mat, j]; 
c = 
            cofactors[mat, j]; 
m . Transpose[c] === 
            IdentityMatrix[Length[m]]*det, {j, 1, Length[mat]}]]checkLaplaceExpansion[mat]

Sparse arrays: "functions" that evaluate by hash table look-up

Q: I have a list of ordered pairs that represent function evaluations. I would
like a quick way to form a "function" f such that for each pair {a,b}, f[a]
will return b.

A: The idea is to form DownValues for f. This uses the internal hashing of Mathematica which will give fast look-up. As an example we will use two lists of ordered pairs of random reals. We do this so we can compare timings.

The particularly neat one-liner to set up the down-values (the hash table) was provided by Allan Hayes.

e1 = Table[{j, Random[]}, {j, (2^12)}]; 
    e2 = Table[{j, Random[]}, {j, (2^16)}]; 
    Timing[(
          f[#1[[1]]] := #1
                [[2]]&)/@e1; ]Timing[
      (g[#1[[1]]] := #1
                [[2]]&)/@e2; ]i1 = Table[Random[Integer, {1, (2^12)}], {1024}];

Now check that accessing a set of fixed size of these down values takes about as long for g as for f.

Timing[Do[Map[f, i1],  {5}]; ]Timing[Do[Map[g, i1],  {5}]; ]

It is actually faster! This mystery will be explained in a moment. A first clue is that if we use an index set spread out more among the possible values for g, then the times are about the same.

\(i2\  = \ Table[Random[Integer, \ {1, 2^16}], \ {1024}]; 

    Timing[Do[Map[g, i2], \ {5}]; ]\)

So we might hypothesize
(i) The internal hash tables are about proportional in length to the number of down values.
(ii) The first index set may have involved elements typically found at the beginning of their hash buckets.
(iii) Typically there are three down values stored per hash bucket because we expect that the well-spread case will need to find on average the middle element, and the first-element case runs twice as fast.

Some measure of confirmation comes from the test below, where we access elements from among the last of the down values. It takes three times as long as accessing from among the earlier-listed values.

\(i3\  = \ Table[Random[Integer, \ {2^16 - 2^12, 2^16}], \ {1024}]; 

    Timing[Do[Map[g, i3], \ {5}]; ]\)

It is clear that the placement of down values within hash buckets must be related to the order in which the down values were created. We can verify this by reversing the order in which we create them and then checking that the timings flip-flop.

\(Map[\((frev[#[\([1]\)]]\  := \ #[\([2]\)])\)&, \ Reverse[e1]]; 

    Map[\((grev[#[\([1]\)]]\  := \ #[\([2]\)])\)&, \ Reverse[e2]]; 


    Timing[Do[Map[frev, i1], \ {5}]; ]\)Timing[Do[Map[grev, i1],  {5}]; ]Timing[Do[Map[grev, i2],  {5}]; ]Timing[Do[Map[grev, i3],  {5}]; ]

In any case, we have constant-access time. One would expect hashing based on monomials rather than integers to be slower by a constant factor.

By the way, timings above show that creation of down values is done in (essentially) O(1) time, and we see from below that it is also roughly constant time to remove them.

ii1  =  Union[i1];
\(\(ii2\  = \ Union[i2];\)
\)
Timing[Map[(f[#] =. ) &, ii1];]
Timing[Map[(g[#] =. ) &, ii2];]

Fast list processing and binary trees

The following problem (slightly altered)  was posed by Richard Gaylord on the news group comp.soft-sys.math.mathematica.

Q: Given a vector of numbers, I wish to form a reordering of this list by taking, at each step, the larger of the two values at the ends of the list, and then moving in one slot from that end to form the new working list. For example, if vec is {18,19,1,11,25,12,22,14} then the result will be {18,19,14,22,12,25,11,1}. How might this be done efficiently?

A: We show two ways to solve this with O(n) complexity, where n is the length of the list. The first method is C-like in that we create a dummy vector and stuff values into it as we find them. This is feasible because we know in advance the length of the result. The second method is Lisp-like in that we build up ever larger lists of the form {old
list, appended value}, and then we flatten it at the end. Timings indicate that we indeed have O(n) behavior.

takeLargestFromEndsC[input_List] := 
      Module[{f = 1, l = Length[input]}, 
        input[[Table[
            If[input[[f]] > 
                input[[l]], f++, 
              l--], {l}]]]]takeLargestFromEndsL[input_List] := 
      Module[{j = 1, k = Length[input], result = {}}, 
        While[j ≤ k, 
          If[input[[j]] > 
              input[[k]], 
            result = {result, 
                input[[j]]}; j++, 
            result = {result, 
                input[[k]]}; k--]; 
          ]; Flatten[result]]

We show some examples.

\(row1K\  = \ Table[Random[Integer, \ 1000], \ {1000}]; 

    row10K\  = \ Table[Random[Integer, \ 1000], \ {10000}]; 


    \(\((tc1K\  = \ takeLargestFromEndsC[row1K])\);\)\  // \ Timing\)(tl1K  =  takeLargestFromEndsL[row1K]);  //  Timing(tc10K  =  takeLargestFromEndsC[row10K]);  //  Timing(tl10K  =  takeLargestFromEndsL[row10K]);  //  Timing

The C-style solution was contributed by Allan Hayes. It runs about twice as fast for large examples, perhaps because it only involves stuffing a pointer to an internal expression in each iteration, whereas the Lisp-like method builds another layer of list as well as inserting a reference, in each iteration. I show the Lisp-like solution primarily to emphasize the tactic of forming nested lists and flattening the result; for large lists this is much faster than use of AppendTo. I should note that deep nesting to the level of tens of thousands of elements will typically crash your version of Mathematica. This has been fixed in our development version.

The previous problem was a warmup exercise. Below is a more challenging variation, posed by Jack Goldberg.

Q: Given n lists (call them l1,...,ln) we wish to select the maximum element from the set of first elements and delete this element from the list it was in. We will repeat this r times. How might we do this efficiently in Mathematica? Further details: at any stage, if there is a tie, we take the entry from the list with the lowest subscript. When a list is empty, the maximum is selected from the remaining lists. Gaylord's problem is a special case in which l1 is the given list and l2  is the given list in reverse order and r == Length[l1].

A: I will show two approaches, one by Michael Trott and one by myself. They show several idiosyncracies of Mathematica programming, hence will get more detailed analysis than most other problems in this Notebook.

My emphasis was on run-time efficiency. To make this reasonable I put in a few assumptions. These are helpful because otherwise we really have special case subproblems e.g. when the number of lists is small, and these can be handled efficiently with code that is not too clever. Moreover different assumptions about distribution of elements might lead to differently optimized algorithms. Our set-up is thus as follows. We have n lists, and we are to build a list of r elements. We assume n is at least, say, 2^6, and also and r>>n. We also assume list elements have uniform random distribution.

Here is a motivation for this problem. It is essentially equivalent to a sort of problem that might be encountered in handling queues: given a bottleneck fed by several pipelines of discrete processes, implement a strategy that iteratively selects and runs the quickest process from among those currently at the front of the pipelines.

The first method uses a fairly simple idea (even I thought of it, I just did not implement it very well). It may be stated as follows.
Keep a list of first elements. Keep a list of indices into the input lists (these record the current front positions, initially all are 1). Iteratively
	Find the largest element in the front positions. Make it the next element in the result under construction.
	Find the first occurrence of that value in the current front positions in the input lists. Increment the index for that list.
	Take the next element from that list, put it into the list of current first elements.

The best implementation I saw of this method was from Michael Trott. Here is his code. I'll explain later why it is so good. Someone on the news group tagged this problem as the Goldberg Variation, hence Michael's naming of the procedure.

\(GoldbergVariationList[lists_, \ n_] := 
\ \ 
      Block[{leadingNumbers, \ leadingPositions, \ listLengths, 

          \ \ \ \ \ \ \ update, \ max, \ pos, \ min}, 
	
        leadingNumbers\  = \ First\ /@\ lists; 
	
        leadingPositions\  = \ Table[1, {Length[lists]}]; 
	
        listLengths\  = \ Length\ /@\ lists; 
	min\  = \ Min[lists]; 
	
        update\  := 
		
          \((max = Max[leadingNumbers]; 
		
            pos = \(Position[leadingNumbers, max, {1}, 1]\)[\([1, 1]\)]; 

            		If[leadingPositions[\([pos]\)] === listLengths[\([pos]\)], 

              			leadingNumbers[\([pos]\)] = min, 
			
              leadingPositions[\([pos]\)] = leadingPositions[\([pos]\)] + 1; 
              
			leadingNumbers[\([pos]\)] = 
                lists[\([pos, leadingPositions[\([pos]\)]]\)]]; 
		max)\); 
        
	Table[update, {n}]
\ \ ]\)

Quick analysis of this method: Adding an element to the result is O(n) to find the largest on the "currently in first position" list. We do this r times. Hence the algorithm is O(r*n). We will verify this later.

The method I show next is different in that we maintain a sorted binary tree of initial elements. With this it is relatively quick to extract the largest or to insert a new element. I think this idea is useful in large part because it shows how one might implement such trees in Mathematica. Though there are pitfalls to this (more on this later).

Here is an outline. First make n pairs of the form {first element, list it lives on}. Put them into a binary tree, sorted by first element (ties will be broken with the second element, enforcing the "take from lowest list" dictum). Iteratively
	Remove pair containing the largest first element. Make it the next element in the result under construction.
	Take next element (if any) from the list it lived on, insert it into tree.

Here is the actual code. I likened the problem to some sort of lottery, hence the choice of name. Before that we define the tree construction and manipulation procedures.

newTree[{}]  :=  {}\(newTree[tlist_List]\  := \ 
      Block[{len, \ len2}, 
	len\  = \ Length[tlist]; 
	
        If\ [len == 1, \ Return\ [{{}, \ tlist[\([1]\)], \ {}}]]; 
	
        len2\  = \ Quotient[len, 2]; 

        	{newTree[Take[tlist, \ len2]], \ tlist[\([len2 + 1]\)], \ 
          newTree[Drop[tlist, \ len2 + 1]]}
	]
\)addToTree[{},  elem_]  :=  { {},  elem,  {} }\(addToTree[tree_, \ elem_]\  /; \ 
        SameQ[tree[\([2, 1]\)], \ elem[\([1]\)]]\  := 
	
      If\ [tree[\([2, 2]\)] < elem[\([2]\)], 

        		{tree[\([1]\)], \ tree[\([2]\)], \ 
          addToTree[tree[\([3]\)], \ elem]}, 

        		{addToTree[tree[\([1]\)], \ elem], \ tree[\([2]\)], \ 
          tree[\([3]\)]}
		]\)\(addToTree[tree_, \ elem_]\  /; \ OrderedQ[{tree[\([2]\)], \ elem}]\  := 
      
	{tree[\([1]\)], \ tree[\([2]\)], \ addToTree[tree[\([3]\)], \ elem]}
        \)\(addToTree[tree_, \ elem_]\  := 
      \ {addToTree[tree[\([1]\)], \ elem], \ tree[\([2]\)], \ tree[\([3]\)]}

    \)removeMaxElement[tree_]  /;  Last[tree]  ==  {}  :=  
      Take[tree, 2]\(removeMaxElement[tree_]\  := \ 
      With[{tmp\  = \ removeMaxElement[tree[\([3]\)]]}, 

        		{{tree[\([1]\)], \ tree[\([2]\)], \ tmp[\([1]\)]}, \ 
          tmp[\([2]\)]}]
\)\(lottery[lists_, \ r_]\  := \ 
      Module[
	{j, \ len = Length[lists], \ firstelems, \ tree, \ elem, \ 
          indices, \ index, \ i = 0}, 
	
        firstelems\  = \ 
          Sort[Table[{First[lists[\([j]\)]], j}, \ {j, \ len}]]; 
	
        tree\  = \ newTree[firstelems]; 
	indices\  = \ Table[1, \ {len}]; 
        
	Table[
		\(i++\); 

          		{tree, \ {elem, j}}\  = \ removeMaxElement[tree]; 
		
          If\ [\((index\  = \(++indices[\([j]\)]\))\) ≤ \ 
              Length[lists[\([j]\)]], 
			
            tree\  = addToTree[tree, \ {lists[\([j, index]\)], j}]]; 
		
          elem, 
		{r}
		]
	]\)

Quick analysis: Initial sorting is O(n*log(n)), each extraction-insertion step is O(log(n)), hence total complexity is O((n+r)*log(n)).

Below are some tests. Among other things they can help to demonstrate the complexity of the algorithms. In all cases I set r to 10000, so doubling this will show the dependence on r at least up to polynomial factors (logarithmic ones can be more subtle to discern).

SeedRandom[1111]test6l = Table[Random[], {(2^6)}, {(2^10)}]; 
    Timing[lot6 = lottery[test6l, 10000]; ]Timing[gvl6 = GoldbergVariationList[test6l, 10000]; ]test9l = Table[Random[], {(2^9)}, {(2^7)}]; 
    Timing[lot9 = lottery[test9l, 10000]; ]Timing[gvl9 = GoldbergVariationList[test9l, 10000]; ]test12l = Table[Random[], {(2^12)}, {(2^4)}]; 
    Timing[lot12 = lottery[test12l, 10000]; ]Timing[gvl12 = GoldbergVariationList[test12l, 10000]; ]

Here are my findings on comparing these algorithms on the examples above and many others not shown, as well as further experimentation with variations of the tree code. In checking on asymptotic behavior, and pondering some timing anomalies in general, I learned several things. First, Michael's code runs like the wind on the first two of the examples shown, in the development version of Mathematica. It is only in the third example that my tree method is demonstrably faster. This is because we now use better internal technology (packed arrays) for storing and accessing elements of arrays of machine numbers. For other types of data the binary-tree approach is superior for much smaller values of n (its behavior is not much dependent on data type so long as it is not packable and order comparisons are relatively fast). Of course it is also the case that prior versions will not be so fast with machine integer data because they do not have the packed array technology. Michael's code also showed an interesting Achilles' heel to array packing. In one pair of examples (not shown), doubling the size of the problem led to a faster result. This was because the smaller problem was below the threshhold of automatic packing, hence used non-packed processing. This has given us cause to re-think some default threshhold settings for use of this packing technology.

The binary-tree method is not amenable to use of this faster technology, so far as I can tell. To implement a binary tree requires alot of nesting and cannot be done with simple tensor structures (e.g. vectors and matrices) that are requisite for the packed-array technology to be utilized.

In the binary-tree algorithm we insert into random positions, but always remove the rightmost element from the binary tree. The tree thus becomes unbalanced over time. We can re-balance the tree if so desired. We simply Flatten it (this will leave it sorted), then Partition into pairs {value, list it came from} and rebuild from this flattened list of pairs. One can avoid the neeed for Partition by using some wrapper other than List, e.g. pair[]. I experimented with this rebalancing and, while it seems to give small improvement in some cases, it more often is disasterous. Hence it was omitted.

Another reason Michael's code is fast is that he avoids the use of subroutines. His update mechanism is defined inside the main routine. To make the tree code run fast I had to avoid use of Module and, to the extent possible, Block constructions as well. These procedural structures are useful when one is doing processing with relatively slow built-in procedures, but they can murder the efficiency of code that is neck-deep in list processing, as, for example, in this problem. Use of DownValue, a standard mode of Mathematica programming, is simply more efficient in this realm. Yet another reasonable alternative is to use, say, Block, but avoid recursion. I implemented a non-recursive version of addToTree, shown below. We give it the HoldFirst attribute to emulate the "call-by-reference" mode for passing the first variable. There is thus now no need to set tree to the result

SetAttributes[addToTree,  HoldFirst]\(addToTree[tree_, \ elem_]\  := \ 
      Block[
	{tmp = tree, \ ll = {}, \ indx}, 
	
        While\ [tmp =!= {}, 
		
          index\  = \ If\ [OrderedQ[{tmp[\([2]\)], \ elem}], \ 3, \ 1]; 
		
          ll\  = \ {ll, index}; 
		tmp\  = \ tmp[\([index]\)]; 
		]; 

        	ll\  = \ Flatten[ll]; 
	
        tree[\([Apply[Sequence, ll]]\)]\  = \ {\ {}, \ elem, \ {}\ }; 
	]
        \)

When I used this I saw no appreciable difference in timing. Hence I did not pursue this with removeMaxElement.

I ran into some further perversity wherein a change of name of a variable local to a Block statement caused a 50% slow-up. This is unfortunate misbehavior, and we are pondering ways to improve on it, but in essence it follows from the "infinite evaluation" mechanism of Mathematica and optimizations related thereto. In other words, it is effectively impossible
to fix completely such problems. In this instance I also found that changing Block to With removed the difficulty, most likely because the latter effectively changes the name of its locally scoped variables.

There is a final aspect to binary tree manipulation that should be mentioned. In the routine above we only add and remove leaves to/from the tree. Whenever we do Part evaluation there is a risk that internal optimizations will fail to detect right away that we need not reevaluate the entire list. This is generally not an issue when all we reevaluate is a bottom-most list containing the altered leaf. More strenuous tree manipulation, with alteration of nodes that are not leaves, might give rise to reevaluation of cascades of nested sublists. (I apologize for my vague attempt at an explanation of when and why lots of reevaluation might take place. In truth I am not fully familiar with the whys and wherefores that give rise to the phenomenon. Suffice it to say that it I have seen worse-than-expected asymptotic complexity of in tree algorithms I wrote that make alterations in the middle of the tree rather than just at the leaves). One way to improve the situation might be to use Hold wrappers, as we did in a list shuffle algorithm earlier in this Notebook. It can be a tricky business to do this with complicated tree manipulation routines. Another way would be to have a built-in expression type that is "inert" to reevaluation, in some sense similar to what we have done with packed tensors. We have experimented with this notion and preliminary results are quite promising.