################################################################### # # # Let f be a rational function in n,k. # # # # Calling sequence: SumRat(f, n, k); # # # # This procedure determines if sum(f, k=0..n-1) is a rational # # function, and if so, it finds that rational function. # # # # For more information, see: # # # # http://arXiv.org/abs/math.CO/0210158 # # # # # # Mark van Hoeij, 2002. # # # ################################################################### # If you find any bugs (which is not unlikely because this code # has not yet been properly tested) please let me know: hoeij@math.fsu.edu # # Apart from possible bugs, the completeness of the algorithm is not # proven, as it depends among other things on an unproven conjecture. # Sum over k = 0 .. n-1 (Note: NOT k=0..n). # SumRat:=proc(f, n, k) # options trace; local F,S,G,s,P,ext; ext := indets([args],{'RootOf','nonreal','radical'}); if ext<>{} then Normalizer:=evala fi; F := convert(f, parfrac, k, ext); # The terms in F. F := [`if`(type(F,`+`),op(F),F)]; (F,P) := group_together(F,n,k); # Would be nice to find more quick tests than just this one: # if member(1, map(nops,F)) then return FAIL fi; S := 0; for G in F do s := SumGroup(G,n,k); # if s = FAIL then return s fi; S := S + s od; S + sum(P,k=0..n-1) end: group_together := proc(F, n, k) # options trace; local Found,Poly,T,D,N,d,m,t,i; Found := { [[k+1, 0, 0, 1, 0, 0]] }; # contains just the rational group Poly := 0; for T in F do if type(T,polynom(anything,k)) then Poly := Poly+T else D := denom(T); N := numer(T); if type(D,`^`) then d := op(1,D); m:=op(2,D) else d := D; m:= 1 fi; i := lcoeff(d,k); if type(i,rational) then # make d monic d := d/i; N := N/i^m fi; t := d, m, N; for i in Found do if same_group(i[1][1],d,k,n,'tr') then Found := (Found minus {i}) union {[op(i),[t, op(tr)]]}; t := 0; break fi od; if t<>0 then Found := Found union { [[t,1,0,0]] } fi fi od; Found := {seq(regroup(i), i=Found)}; Found, Poly end: # Put each multiplicity seperately. regroup := proc(G) local T,i,j; T := {seq(i[2],i=G)} minus {0}; # multiplicities seq([`if`(G[1][2]=0, subsop(2=j,G[1]), NULL), seq(`if`(i[2]=j,i,NULL),i=G)], j=T) end: # Some simple properties that should be the same for two polynomials # in the same group. properties := proc(f, k, n) [degree(f, k), degree(l, n), degree(f,{k,n}), degree(lcoeff(f,k),n), indets(f,{'name','RootOf','radical'}) minus {n}] end: # Take coefficients w.r.t. everything but vars cfs := proc(F, vars) local v,l,j; v := indets(F,{'RootOf','radical'}); if v = {} then {seq(op([collect(i,vars,Re),collect(i,vars,Im)]), i = [coeffs(F, indets(F,'name') minus vars)])} else # RootOfs not contained in other RootOfs v := v minus indets( [seq(op(1,j),j=v)] ,{'RootOf','radical'}); l := proc() local a; a end: procname(subs({seq(j=l(),j=v)},F), vars) fi end: # Let p(n) be an algebraic function in n. # Then k - p(n) is in the same group as # r1*k + r2*1 + r3*n - p(n) where r1,r2,r3 are rational numbers. # same_group := proc(f,g,k,n,tr) # options trace; local r1,r2,r3,c,c2,zero,vars,S,B; if properties(f,k,n) = properties(g,k,n) then c := Normalizer(lcoeff(f,k)/lcoeff(g,k)); if not has(c,n) then zero := subs(k=r1*k + r2*1 + r3*n, f) - c * r1^degree(f,k) * g; vars := {r1,r2,r3}; zero := cfs(evala(Expand( zero )), vars); # if {degree(f,k),degree(f,{k,n})}={1} and indets(f, # {'name','RootOf','radical','nonreal'}) minus {k,n} = {} then if f = k+1 then # Group = k. zero := zero union {r1-1} fi; S := [solve(zero, vars)]; c := [r1,r2,r3]; S := [seq(`if`(type(subs(i,c),list(rational)),i,NULL),i=S)]; if S<>[] then c := subs(S[1], c); B := evalb( type(c, list(rational)) ); if B and nops(S)=2 then if evala(subs(S[1],r1)+subs(S[2],r1))<>0 then error "should not happen" else c2 := subs(S[2], [r1,r2,r3]); if not type([op(c),op(c2)],list(rational)) then return false fi; if c[3] + min(0,c[1]) < c2[3] + min(0,c2[1]) then c, c2 := c2, c fi; if nargs<7 then procname(f,f,k,n,'S',0,0); if S[1]=-1 then c := c2 elif S[1]<>1 then error "should not happen" fi fi fi fi; tr := c; return true fi fi fi; false end: NormalizeTerm := proc(i, n, k, RAT) # options trace; local c, j, S, T1, T2,s; c := `if`(RAT, ceil(i[5]), floor(i[5]/i[4])); if c = 0 then # already normalized [op(i), 0] # 0 = impact of normalization on the sum else j := subsop(5 = i[5] - c*i[4], i); T1 := i[3]/i[1]^i[2]; T2 := subs(k=k+n, T1); if c > 0 then S := add(subs(k=s, T2-T1), s=-c..-1) else S := add(subs(k=s, T1-T2), s=0..-c-1) fi; subs(k=k-c,[op(j), S]) fi end: CombineTerms := proc(L, n, k) # options trace; local i,j,h,a,b; for i to nops(L) do a := L[i]; if Middle_in_Symgroup(a, n, k) then if not `dsolve/diffeq/higherorder/iszero`(a[3]-subs(k=n-1-k,a[3])) then b := collect((a[3] + subs(k=n-1-k,a[3]))/2,{k,n},Normalizer); return procname(subsop(i = `if`(b=0,NULL,subsop(3=b,a)), L), n, k) fi; for j to nops(L) do b := L[j]; if j<>i and a[4]=b[4] and a[6]=b[6] and (a[5]-b[5])/a[4]=1/2 then for h to nops(L) do if L[h][4] = 1/2 * a[4] then return TreatMiddle(a,[seq(`if`(h=i,NULL,L[h]),h=1..nops(L))], n, k) fi od fi od fi; for j from i+1 to nops(L) do b := L[j]; if a[4..6] = b[4..6] then a := [a[1], a[2], a[3] + b[3] * (lcoeff(a[1],k)/lcoeff(b[1],k))^a[2], a[4], a[5], a[6]]; if testeq(a[3]) and Normalizer(a[3])=0 then a := NULL fi; return procname(subsop(j=NULL, i=a, L), n, k) fi od od; L end: Middle_in_Symgroup:=proc(a, n, k) # options trace; evalb( irem(degree(a[1],k),2) = 0 and `dsolve/diffeq/higherorder/iszero`(a[1]-subs(k=n-1-k,a[1])) ) end: TreatMiddle := proc(a,Rest,n,k) local a3,t; global `SumRat/side`; # options trace; a3 := (a[3]+subs(k=n-1-k,a[3]))/2; t:=a3/a[1]^a[2]; `SumRat/side`:=`SumRat/side`-subs(k=-1/2,t)+subs(k=(n-1)/2,t); CombineTerms([ [op(subs(k=(k-1)/2,[a[1],a[2],2*a3])), a[4]/2,a[5]-a[4]/2,a[6]], [op(subs(k=k-1/2, [a[1],a[2], -a3])), a[4],a[5]-a[4]/2,a[6]] ,op(Rest) ], n, k) end: SumGroup := proc(G, n, k) # options trace; local L, i, S, RAT; global `SumRat/side`; RAT := evalb( G[1][1] = k+1 ); if RAT then # Make r3 = i[6] at least -1/2. L := [seq(`if`(i[6] >= -1/2, i, [subs(k = n-1-k, -i[1]), i[2], -i[3], 1, -1-i[5], -1-i[6]]) ,i = G[2..-1])] else # Make the speed r1 = i[4] positive: L := [seq( `if`(i[4]>0, i, [op(evala(Expand(subs(k = n-1-k, i[1..3])))), -i[4], i[5]-i[4], i[6]+i[4]] ), i = G)]; fi; # Now normalize r2 = i[5], i.e. move it in [0, r1). L := map(NormalizeTerm, L, n, k, RAT); S := add(i[7], i=L); L := [seq(evala(Expand(i[1..6])),i=L)]; `SumRat/side` := 0; L := CombineTerms(L, n, k); S := normal(S+`SumRat/side`); if ZeroSum(L, n, k) then S else S + Sum(add(i[3]/i[1]^i[2], i=L), k=0..n-1) fi end: ZeroSum := proc(L, n, k) # options trace; # Quick test. if nops(L) = 1 or not CheckSmallN(L, n, k) then return false fi; # Should do additional checking here, not yet implemented. # It is possible to prove zero-sum by checking only finitely # many k. true end: CheckSmallN := proc(L, n, k) local T,i; T := add(i[3]/i[1]^i[2], i=L); `dsolve/diffeq/higherorder/iszero`( [seq(add(subs(n=i,k=s,T),s=0..i-1), i=1..4)]) end: # t:=1/(k*(2*n-k)+1); # F5 := t + subs(k=2*n-k,t) - subs(k=2*k,t) - subs(k=2*k-1,t); # F:=normal(F5); # F := (4*k-2*n-1-2*n*k+16*n^2*k^2-40*n*k^3+k^2+24*k^4+38*k^2*n-12*n^2*k-28*k^3) # /(-2*n*k+k^2-1)/(-4*n*k+4*k^2-1)/(-2*n*k+2*k^2-2*k+n); # F := (k-((n-1)/2))/( (k-((n-1)/2))^2 + n); # recently fixed bug. # SumRat(F,n,k);