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] ====