# The code below proves Worpitzky's identity by induction.
#
#
#
#
#
interface(screenwidth=68):
eulerian := proc(n,k)
 if n=k or type(k,numeric) and k<0 then 0;
 elif n=1 and k=0 then 1;
 else 'eulerian'(n,k);
 fi;
end:

`expand/eulerian` := proc(n,k)
 if type(n,`+`) then
  (k+1)*eulerian(n-1,k) + (n-k)*eulerian(n-1,k-1)
 else eulerian(n,k)
 fi;
end:

distribute := proc(a)
 if type(a,`*`) and type(op(1,a),`+`) then
  map(proc(a,b) a*b end,op(1,a),a/op(1,a))
 else
  a
 fi;
end:

`expand/Sum` := proc(s) local a,b;
 if type(s,Sum(`+`,`=`)) then
  a := op(1,s); b := op(2,s);
  map(Sum,a,b)
 else s
 fi
end:

remove_constant := proc(s) local a,b,v,av;
 if type(s,Sum(`*`,`=`)) then
  a := op(1,s); b := op(2,s); v := op(1,b);
  av := map(proc(a,v) if has(a,v) then a else 1 fi end,a,v);
  a/av*Sum(av,b);
 else
  s;
 fi;
end:
 
 
formula := x^n=Sum(eulerian(n,k)*binomial(x+k,n),k=0..n-1):
lprint(`Prove  Worpitzky's identity.`);
lprint(` Prove`);
print(formula);
lprint(` for n>0 .`);
lprint(` Check that the case n=1 is true.`);
L1 := eval(subs([Sum=sum,n=1],formula));
lprint(` By the inductive assumption, assume the formula`);
lprint(` is true for n=n and try to prove it for n=n+1 .`);
Assume := rhs(formula)=x^n;
lprint(` Evaluate the right hand side at n=n+1 .`);
L2 := subs(n=n+1,rhs(formula));
lprint(` Use the eulerian number recurrance relation.`);
L3 := subs(eulerian(n+1,k)=expand(eulerian(n+1,k)),L2);
lprint(` Distribute the binomial over the sum.`);
L4 := map(distribute,L3);
lprint(` Expand the sum into 2 sums.`);
L5 := `expand/Sum`(L4);
lprint(` Remove the last term of the first sum since eulerian(n,n)=0 .`);
L6 := subsop(1=subs(0..n=0..n-1,op(1,L5)),L5);
lprint(` Remove the first term of the second sum since eulerian(n,-1)=0 .`);
L7 := subsop(2=subs(0..n=1..n,op(2,L6)),L6);
lprint(` Let k=j+1, in the second sum.`);
L8 := subsop(2=Sum(subs(k=j+1,op(1,op(2,L7))),j=0..n-1),L7);
lprint(` Change j to k.`);
L9 := subs(j=k,L8);
lprint(` Combine the sums.`);
L10 := Sum(op(1,op(1,L9))+op(1,op(2,L9)),op(2,op(1,L9)));
lprint(` Expand the binomial coefficients.`);
L11 := subs([binomial(x+k+1,n+1)=expand(binomial(x+k+1,n+1)),
            binomial(x+k  ,n+1)=expand(binomial(x+k  ,n+1))],L10);
lprint(` Normalise inside the sum.`);
L12 := normal(L11);
lprint(` Move x outside the sum.`);
L13 := remove_constant(L12);
lprint(` Use our inductive assumption.`);
L14 := subs(Assume,L13);
lprint(` Combine the powers of x.`);
L15 := simplify(L14,power);
lprint(` Check that this is equal to the    `);
lprint(` left hand side evaluated at n=n+1 .`);
evalb(L15=subs(n=n+1,x^n));


### For use by authors only
### Experimental editor for the Maple Form Interface
###


### For use by authors only
### Experimental editor for the Maple Form Interface
###


### For use by authors only
### Experimental editor for the Maple Form Interface
###


### For use by authors only
### Experimental editor for the Maple Form Interface
###