On 3/14/07, Kostas Oikonomou <ko at research.att.com> wrote:
> Stavro,
>
> If so, could you provide this trivial implement ion
> KostasSum :-) -> nested sum using macros?
Kosta, try this:
expandsum(expr,[ranges])::=expandsum_helper(expr,ranges)$
/* macro definition -- main entry, but not useful for recursion */
expandsum_helper(expr,ranges):=
if ranges=[] then expr
else block([range: first(ranges)],
if length(range)#3 then
error("Index range must be triple:",
"[varname, lowval, highval]",
range),
funmake(sum,[expandsum_helper(expr,rest(ranges)),
range[1],
range[2],
range[3] ]))$
Test cases:
expandsum(23) => 23
expandsum(f(i),[i,0,0]) => f(0)
expandsum(f(i),[i,1,0]) => 0
expandsum(f(i,j),[i,1,j],[j,1,n]) => sum(sum(f(i,j),j,1,n),i,1,j)
(probably not what you want)
expandsum(f(i,j),[j,1,n],[i,1,j]) => sum(sum(f(i,j),i,1,j),j,1,n)
(if you want the nesting to go the other way, change
expandsum to
... ::=expandsum_helper(expr,reverse(ranges)) )
expandsum(i+j,[j,1,n],[i,1,j]) => sum(sum(j+i,i,1,j),j,1,n)
expandsum(i+j,[j,1,n],[i,1,j]), simpsum
=> (2*n^3+3*n^2+n)/4+(n^2+n)/4
expandsum(i*j,[j,1,n],[i,1,j]) => sum((sum(i,i,1,j))*j,j,1,n)
(note simplification)
expandsum(i*j,[j,1,n],[i,1,j]), simpsum
=> ((n^4+2*n^3+n^2)/4+(2*n^3+3*n^2+n)/6)/2
factor(%) => n*(n+1)*(n+2)*(3*n+1)/24
Now let's use the Gosper sum algorithm...
substnusum(expr):=ev(subst('nusum,nounify('sum),expr),eval)$
/* This should be made more bulletproof, so that it sums what it
can and leaves the rest alone, rather than giving an error... */
expandsum(j!,[j,1,n],[i,1,j]) => sum(j*j!,j,1,n)
substnusum(%) => (n+1)!-1
All the work above is being done by the existing "sum" code.
Now let's write the inverse:
/* Converts all sums (even non-nested ones) to expandsum form */
contractsums(ex):=
block([summ:nounify('sum)],
if atom(ex) then ex
elseif inpart(ex,0)=summ then
block([build:[]],
while inpart(ex,0)=summ do
(build:cons([inpart(ex,2),inpart(ex,3),inpart(ex,4)],build),
ex:inpart(ex,1)),
funmake('expandsum,cons(ex,reverse(build))))
else map(contractsums,ex))$
expandsum(f(i,j,k),[k,1,n],[j,1,k],[i,1,j]) =>
sum(sum(sum(f(i,j,k),i,1,j),j,1,k),k,1,n)
contractsums(%) =>
expandsum(f(i,j,k),[k,1,n],[j,1,k],[i,1,j])
ev(%) => sum(sum(sum(f(i,j,k),i,1,j),j,1,k),k,1,n)
Hope this helps,
-s