more Goldberg Variations

*To*: mathgroup at smc.vnet.net*Subject*: [mg4064] more Goldberg Variations*From*: Vandemoortele CC Group R&D Center <wouter.meeussen.vdmcc at vandemoortele.be>*Date*: Thu, 30 May 1996 02:50:09 -0400*Sender*: owner-wri-mathgroup at wolfram.com

I tried to post aan amendment to my previous note on this subject, but did not seem to get onto the bulletin board. M. Berth seems to have had similar problems, and requested me to post his contribution for him. Yesterday, I got a clarification from Allan Hayes [mg4015] who did get on the bulletin board. Since Berth made no mention of this, I will also reproduce it here. Allan Hayes wrote: start--of--quote > This is also done by the function Compositions from the standard > package DiscreteMath`Combinatorica` : > > <<DiscreteMath`Combinatorica` > > Compositions[3,3] > {{0, 0, 3}, {0, 1, 2}, {0, 2, 1}, {0, 3, 0}, {1, 0, 2}, {1, 1, 1}, > {1, 2, 0}, {2, 0, 1}, {2, 1, 0}, {3, 0, 0}} > >Allan Hayes >hay at haystack.demon.co.uk > end--of--quote Now follows what M. Berth asked me to post for him: START-OF-QUOTE I have come across this problem before (in the context of mixed partial derivatives of a function), it can be stated as: Generate a list of all multiindices J=(i1, i2, ..., ip), 0 <= ik <= s such that i1 + i2 + ... + ip = s. Let's call p the length and s the degree of J. Then it is easy to find the multiindices of length p and degree 1 -- it's just the "unit vectors" (1,0,...,0), (0,1,0,...,0) etc. The multiindices of length p and degree 2 are then found by adding unit vectors. So for each multiindex of degree k we get several multiindices of degree k+1 by adding suitable unit vectors. This gives a recursive definition, like in the following program: (** Multiindices **) (* Unit vectors are the primitive multiindices *) units[n_]:=IdentityMatrix[n] (* lnze (read "last non-zero element") determines the index of the rightmost nonzero element *) zeroQ[x_]:=(x==0) lnze[{x___,y___?zeroQ}]:= Length[{x}] (* next[J] generates the multiindices J,i with degree |J|+1 *) next[{x__?zeroQ}] := units[Length[{x}]] next[l_List] := Map[Plus[#,l]&, Take[units[Length[l]],{lnze[l],Length[l]}]] (* MultiIndices[s,p] gives a list of all multiindices of degree s and length p *) MultiIndices[1,p_] := units[p] MultiIndices[s_,p_]:= Flatten[Map[next,MultiIndices[s-1,p]],1] /; (s>1) && (p>=1) (** The End **) So... this works, but how do you choose these "suitable" unit vectors, such that no multiindex is created twice? The function next takes a multiindex J of degree s and gives the list of multiindices of degree s+1 generated from J by adding those unit vectors. If the rightmost index is greater than zero, next will just add (0,0,...,1). If there are one or more zeros in the right part of J (like (..., r, 0, ..., 0)), then next generates a list of several multiindices containing (..., r+1, 0,...,0), (...,r,1,0,...,0), (...,r,0,1,0,...,0) , ... , (...,r,0,...,0,1). For multiindices of length 2, one can show the effect of next[] by the following picture: (0,0) (order 0) (1,0) (0,1) (order 1) (2,0) (1,1) (0,2) (order 2) (3,0) (2,1) (1,2) (0,3) (order 3) ... (It looks a bit like the Pascal triangle of binomial coefficients... ) Each entry can be generated by adding a unit vector to either its upper or its upper left neighbour. So next[] generates the lower right neighbour(s). It is now possible to give an inductive proof for the correctness of our approach --- begin sketch of proof --- We prove the correctness for arbitrary length p by induction over the order s. For order s=1 the program gives the correct answer -- the unit vectors. Suppose, the answer is correct for an arbitrary order s, i.e. MultiIndices[s,p] gives all multiindices of order s and length p. Then, for each multiindex K of order s+1 we can find a _unique_ predecessor J among the multiindices of order s by subtracting 1 from the rightmost nonzero element of K. The function next[] takes one such predecessor J and generates all successors K1, K2, ... None of the multiindices K of order s+1 will be generated twice, because each K can be generated only from its unique predecessor. --- end sketch of proof --- I hope this helps to explain the program a bit. -- Matthias Matthias Berth Institute of Mathematics and Computer Science University of Greifswald, Germany berth at rz.uni-greifswald.de END-OF-QUOTE my own second trial to generate a flexible iterator now follows: A generator for a set of k indices as iterators for Table or Sum, subject to the condition that their sum be equal to n can be realised by : index[k_,n_]:=Module[{t}, Sequence@@ ReplacePart[t=Transpose[ {Array[a,k], Table[0,{k}], FoldList[#-#2&,n,Array[a,k-1]]} ], t[[k,3]],{k,2}] ] it produces a sequence of iterators : in:= index[3,5] out:=Sequence[{a[1], 0, 5}, {a[2], 0, 5 - a[1]}, {a[3], 5 - a[1] - a[2], 5 - a[1] - a[2]}] For some unclear reason, this works only properly inside a 'Table' or 'Sum' when an 'Evaluate' is wrapped around the 'index' function. We end up with what I think Goldberg wanted: in[]:= Sum[Times@@({x,y,z}^Array[a,3]) Multinomial[Sequence@@Array[a,3]], Evaluate[index[3,5]]]//InputForm out[]:= x^5 + 5*x^4*y + 10*x^3*y^2 + 10*x^2*y^3 + 5*x*y^4 + y^5 + 5*x^4*z + 20*x^3*y*z + 30*x^2*y^2*z + 20*x*y^3*z + 5*y^4*z + 10*x^3*z^2 + 30*x^2*y*z^2 + 30*x*y^2*z^2 + 10*y^3*z^2 + 10*x^2*z^3 + 20*x*y*z^3 + 10*y^2*z^3 + 5*x*z^4 + 5*y*z^4 + z^5 in[]:= %//Simplify out[]:= (x + y + z)^5 Too bad that the implicit variables a[..] need definition inside the index[_,_] function, so that we are not at liberty to choose a name for the iterators. If a iterator name could be passed, it would even be better.( Of course, the iterators need a name: they should be 'addressable' from the body of the Sum or Table function.) I checked the speed of Table[Array[a,7],Evaluate[ index[7,7]]~Flatten~6 versus Compositions[7,7] and concluded they both run at about the same speed (7 sec) and give identical results. I hope this gets through the 'daemons' barrior, Wouter. ==== [MESSAGE SEPARATOR] ====