# Recc.ml or Math/elim/solv0.m 1993/12/10,14 ################## try the examples below ### Example 1 (Apery's recurrence). getRecurrence(2,[-n,-n,1+n,1+n,1,1,1],n); ### The output [[1,0,-1],[c1,c2,c3],[d1,d2,d3]] means the given ### hypergeometric sum $f(n):={}_4 F_3(-n,-n,1+n,1+n;1,1,1;1)$ satisfies ### the recurrence ### c1*d1*f(n+1) + c2*d2*f(n) + c3*d3*f(n-1)=0. ### ### Example 2. getRecurrence(4,[a1+n,a2,a3,a4,a5,b1,b2,b3,b4],n); ### You will get a 4th order difference equation with respect to n ### satisfied by the hypergeometric sum ### $ {}_5 F_4 (a1+n,a2,a3,,a4,a5;b1,b2,b3,b4;1)$ ################## Notice that ## x, Dx, g[i], ## a[i],b[i],t (for hg), Verbose, L ## are the gloval variables. ###################################################### with(linalg): readlib(factors): ### Messages. print(); print(`Loading getRecurrence(rank,[parameters],var) : Version 1993-12-16.`); print(`The prototype of this package was written by the computer algebra system kan (133.30.90.19 public/Kan). This is the maple version.`); print(); Verbose:=false: ### Put true if you like verbose messages. #Verbose:=true: ### Put false if you don't like verbose messages. mycoeff:=proc(f,k) ## ex. mycoeff(n*Dx + 3,1) --> n local ff; ff := expand(f); if k = 0 then RETURN(subs(Dx=0,ff)); else RETURN(coeff(ff,Dx,k)); fi; end : myLength:=proc(f) ## ex. myLength([1,2,5]) --> 3 local i,s; s := 0; for i in f do s := s+1; od; RETURN(s); end : getRule0:=proc(f) local deg,topc,ans,i; deg := degree(f,Dx) ; topc := -factor(mycoeff(f,deg)) ; ans := 0 ; for i from (deg-1) by -1 to 0 do ans := ans + mycoeff(f,i)*g[i] ; od ; ans := ans / topc ; RETURN([deg,ans]); end : getRule:=proc(b) # getRule(bbb); give the output of sequence local rule,n ,tmp; n := myLength(b) ; if Verbose then print(`We are in the getRule.`,n); fi: tmp := getRule0(b[1]) ; # tmp = [i, f] rule := {g[tmp[1]] = numer(tmp[2])/factor(denom(tmp[2]))} ; for i from 2 by 1 to n do tmp := getRule0(b[i]) ; tmp := subs(rule,tmp) ; rule := rule union {g[tmp[1]] = numer(tmp[2])/factor(denom(tmp[2]))} ; ##print(rule); if Verbose then print(i/n); fi: od; RETURN(rule); end: tog:=proc(f) ## Dx^2 + n Dx + 1 ---> g[2] + n g[1] + g[0] local deg,topc,ans,i; deg := degree(f,Dx) ; ans := 0 ; for i from deg by -1 to 0 do ans := ans + mycoeff(f,i)*g[i] ; od ; RETURN(ans) ; end : getMat:=proc(m,n,aaa,rule) ## g[0], ..., g[m]: ## aaa is the output of compose 0 ## n is the rank of the difference eq. n and m should be equal. ## c[i] is global ## example read `data1.ml`; getMat(2,2,aaa,getRule(bbb)); linsolve(",[0,0,0]); ## [(n+1)^6, -(2 n + 3) (17 n^2 + 51 n + 39) , 1] local i,vec,tmp,ans,j,gi,ppp,fff; if Verbose then print(`We are computing the matrix.`,m+1,n+1); fi: ans := [ ] ; tmp := 0; for i from 0 by 1 to n do fff :=subs(rule,tog(aaa[i+1][3])); tmp := tmp + expand(c[i]*fff); od: #checkDegree(tmp,m); for i from 0 by 1 to m do vec := [ ] ; gi := coeff(tmp,g[i],1); for j from 0 by 1 to n do vec := [op(vec),factor(coeff(gi,c[j],1))]; od : ans := [op(ans),vec] ; od : ## ans is -------- n+1 ---------- ## I ## m+1 ## I ## -------- n+1 ---------- matrix RETURN(ans); end : # for debug. check the g[i] checkDegree:=proc(f,m) local i; if Verbose then print(`Checking the tmp in getMat`); fi: for i from m+1 by 1 to 60 do if not (degree(f,g[i]) = 0) then print(`error `,i); fi: od: if Verbose then print(`done`); fi: end: factorvecs:=proc(mat) local i,j,ggcd,ggcd2,gcdvec,gcdvec2,m,n; ## mat is -------- n ---------- ## I ## m ## I ## -------- n ---------- matrix m:=myLength(mat); n:=myLength(mat[1]); gcdvec:=[]; for j from 1 by 1 to n do ggcd:=0; ggcd2:=0; for i from 1 by 1 to m do ggcd:=gcd(ggcd,numer(mat[i][j])); ggcd2:=gcd(ggcd2,denom(mat[i][j])); od: if ggcd = 0 then ggcd:=1 fi; if ggcd2 = 0 then ggcd2:=1 fi; gcdvec:=[op(gcdvec),factor(ggcd/ggcd2)]; od: gcdvec2:=[]; for i from 1 by 1 to m do ggcd:=0; ggcd2:=0; for j from 1 by 1 to n do ggcd:=gcd(ggcd,numer(mat[i][j])); ggcd2:=gcd(ggcd2,denom(mat[i][j])); od: if ggcd = 0 then ggcd:=1 fi; if ggcd2 = 0 then ggcd2:=1 fi; gcdvec2:=[op(gcdvec2),factor(ggcd/ggcd2)]; od: factorVecs:=[gcdvec,gcdvec2]; RETURN([gcdvec,gcdvec2]); end: factorout:=proc(mat) local i,j,gcdvec,gcdvec2,m,n,tmp,row,ans,ans2; m:=myLength(mat); n:=myLength(mat[1]); tmp:=factorvecs(mat); gcdvec:=tmp[1]; ans:=[]; for i from 1 by 1 to m do row := [ ]; for j from 1 by 1 to n do row:=[op(row), factor(mat[i][j]/gcdvec[j])]; od: ans := [op(ans),row]; od: gcdvec2:=gcdvec; tmp:=factorvecs(ans); gcdvec:=tmp[2]; ans2:=[]; for i from 1 by 1 to m do row := [ ]; for j from 1 by 1 to n do row:=[op(row), factor(ans[i][j]/gcdvec[i])]; od: ans2 := [op(ans2),row]; od: RETURN([ans2,gcdvec2]); end: ################## functions for multiplication myDiffDx:=proc(f,k) if k=0 then RETURN(f) else RETURN(diff(f,Dx$k)); fi end: myDiffx:=proc(f,k) if k=0 then RETURN(f) else RETURN(diff(f,x$k)); fi end: myfac0:=proc(n,k) local r,i; if k = 0 then RETURN(1); else r:=n; for i from 1 by 1 to (k-1) do r:=r*(n-i); od: RETURN(r); fi end: ppMult:=proc(f,g) local n, r,k,q,p; q:=expand(f) ; p:=expand(g); n := degree(q,Dx); r := 0; for k from 0 by 1 to n do r := r + (1/myfac0(k,k))*myDiffDx(q,k)*myDiffx(p,k) ; od : RETURN(expand(r)); end: ppMult00:=proc(f,g) local n, r,k,q,p; q:=expand(f) ; p:=expand(g); n := degree(q,Dx); r := 0; for k from 0 by 1 to n do r := r + (1/myfac0(k,k))*myDiffDx(q,k)*subs(x=0,myDiffx(p,k)); od : RETURN(expand(r)); end: bench1:=proc(n) local r,i; r:=1; for i from 1 to n do r:= ppMult(x*Dx+2,r); od: RETURN(r); end: #### to get tables compose0:=proc(param,stepup,m) ## compose0(n,[1,(n+1),(Dx+1)],3) ## only for up local ans0,i,r,step; if Verbose then print(`we are in compose0`); fi: step:=stepup[3]; ans0:=[[1,1,1]]; r:=1; for i from 1 by 1 to m do r:=ppMult00(subs(param=param+1,r),step); ans0:=[op(ans0),[factor(subs(param=param+1,ans0[i][1]))*stepup[1], factor(subs(param=param+1,ans0[i][2]))*stepup[2], r]] ; if Verbose then print(i/m); fi: od: RETURN(ans0); end: compose0down:=proc(param,stepdown,m) ## compose0down(n,[1,(n+1),(Dx+1)],3) ## only for up local ans0,i,r,step,ra,rb; if Verbose then print(`we are in compose0down`); fi: step := stepdown[3]; ans0:=[]; r:=1; ra:=1; rb:=1; for i from 1 by 1 to m do r:=ppMult00(subs(param=param-1,r),step); ra:=factor(subs(param=param-1,ra))*stepdown[1]; rb:=factor(subs(param=param-1,rb))*stepdown[2]; ans0:=[op(ans0),[ra,rb,r]] ; if Verbose then print(i/m); fi: od: RETURN(ans0); end: sequence:=proc(ell,m) local ans1,i,dd; if Verbose then print(`we are in sequence`); fi: ans1:=[subs(x=0,ell)]; dd:=degree(ans1[1],Dx); ## dd+i <= m for i from 1 by 1 to (m-dd) do ans1:=[op(ans1),subs(x=0,ppMult(Dx^i,ell))]; od: RETURN(ans1); end: #### to get contiguity deltaToD:=proc(g,t) local n,i,result,f; f:=expand(g); n:=degree(f,t); result:=subs(t=0,f); for i from 1 by 1 to n do result:=result+coeff(f,t,i)*pPower(x*Dx,i); od: RETURN(result); end: pPower:=proc(f,p) local r,i; if p=0 then RETURN(1) fi; r:=1; for i from 1 by 1 to p do r:=ppMult(r,f); od: RETURN(r); end: hg1:=proc(n) local r,i; r:=t; for i from 1 by 1 to n do r:=r*(t+b[i]-1); od: RETURN(r); end: hg2:=proc(n) local r,i; r:=1; for i from 1 by 1 to n do r:=r*(t+a[i]); od: RETURN(r); end: upa:=proc(k) RETURN([1,a[k],t+a[k]]); end: downa:=proc(k,p) #for p F p-1 local ss,c2,tt; ss :=factor(hg2(p)/(t+a[k])) ; c2 :=subs(t= -a[k],hg1(p-1)); tt :=factor((hg1(p-1)-c2)/(t+a[k]));; ss :=factor(subs(a[k] = a[k]-1,ss)); tt :=factor(subs(a[k] = a[k]-1,tt)); c2 :=factor(subs(a[k] = a[k]-1,c2)); c2 := c2/(a[k]-1); RETURN([1,-c2,tt-x*ss]); end: downb:=proc(k) RETURN([1,b[k]-1,t+b[k]-1]); end: upb:=proc(k,p) local ss,c2,tt; tt :=factor(hg1(p-1)/(t+b[k]-1)) ; c2 :=subs(t= -b[k]+1,hg2(p)); ss :=factor((hg2(p)-c2)/(t+b[k]-1));; ss :=factor(subs(b[k] = b[k]+1,ss)); tt :=factor(subs(b[k] = b[k]+1,tt)); c2 :=factor(subs(b[k] = b[k]+1,c2)); RETURN([b[k],c2,x^(-1)*tt-ss]); end: ############################################ getUp:=proc(p,i) # p F p-1, if i<= p then RETURN(upa(i)); else RETURN(upb(i-p,p)) fi; end: getDown:=proc(p,i) if i<=p then RETURN(downa(i,p)); else RETURN(downb(i-p,p)) fi; end: getSubs:=proc(rule,p) # p F p-1 [1,2,3,n,m] -->{a[1]=1,a[2]=2,a[3]=3, ## b[1]=n,b[2]=m} local i,ans; ans := {}; for i from 1 by 1 to p do ans := ans union {a[i]=rule[i]}; od: for i from 1 by 1 to (p-1) do ans := ans union {b[i] = rule[i+p]}; od: RETURN(ans); end: getL1:=proc(param,rule,p,mysign) ## getL1(n,[-n,n+1,1],2,+1)---> [a,b,ell]; p F p-1 # param = param + mysign local ruleup,n,sa,wa,i,co,crule,j,pp,pp0,csubs,L1a,L1b,L1; ruleup:=subs(param=param+mysign,rule); n:=myLength(rule); sa:=[]; wa:=1; for i from 1 by 1 to n do sa:=[op(sa),expand(ruleup[i]-rule[i])]; wa:=wa + abs(sa[i]); od: if Verbose then print(sa); fi: co:=1; crule:=rule; L1:=1; L1a:=1; L1b:=1; for i from 1 by 1 to n do if sa[i] > 0 then for j from 1 by 1 to sa[i] do pp0 := getUp(p,i); csubs := getSubs(crule,p); L1a := L1a*factor(subs(csubs,pp0[1])); L1b := L1b*factor(subs(csubs,pp0[2])); pp:=expand(deltaToD(pp0[3],t)); pp:=subs(csubs,pp); L1:=ppMult(pp,L1); if Verbose then print(co/wa); fi: co :=co+1; crule:=subsop(i=crule[i]+1,crule); od; fi: if sa[i] < 0 then for j from -1 by -1 to sa[i] do pp0 := getDown(p,i); csubs := getSubs(crule,p); L1a := L1a*factor(subs(csubs,pp0[1])); L1b := L1b*factor(subs(csubs,pp0[2])); pp:=expand(deltaToD(pp0[3],t)); pp:=subs(csubs,pp); L1:=ppMult(pp,L1); if Verbose then print(co/wa); fi: co :=co+1; crule:=subsop(i=crule[i]-1,crule); od; fi: od: RETURN([L1a,L1b,L1]); end: zeroVec:=proc(size) local ans,i; ans:=[]; for i from 1 by 1 to size do ans:=[op(ans),0]; od: RETURN(ans); end: myReverse:=proc(list) local n,ans,i; n:=myLength(list); ans:=[]; for i from n by -1 to 1 do ans:=[op(ans),list[i]]; od: RETURN(ans); end: getDifferenceOp:=proc(rup,rdown) local i,ans; ans:=[]; for i from rup by -1 to 0 do ans:=[op(ans),i]; od: for i from -1 by -1 to -rdown do ans:=[op(ans),i]; od: RETURN(ans); end: getBestDecomposition:=proc(r,degL1,degL1down) local dd,i,tmp; dd:=r*degL1down; choice := [0,r]; for i from 0 by 1 to r do tmp:=max(i*degL1, (r-i)*degL1down); if tmp < dd then choice:=[i,r-i]; dd:=tmp; fi; od: RETURN(choice); end: invalid:=proc(clist) local size,i; size:=myLength(clist); for i from 1 by 1 to size do if (clist[i][1] = 0) or (clist[i][2] = 0) then RETURN(true); fi: od: RETURN(false); end: getRecurrence:=proc(r,rule,param) ## r is an expected rank of difference eq. ex. 2 ## rule ex. [-n,-n,n+1,n+1,1,1,1]; ## param ex. n ## They should be global variables when you want to call me manually. # local L,L1,L1down,mat,clist,alist,subsrule,factorVec ## true local variables local degL1,degL1down,rup,rdown,ans,base,i,exx,tmp,p; p := (myLength(rule)+1)/2; base := p-2; ## g[0],...,g[base] is the base. cf.getRule(sequence(L,10) ## p p F p-1 ex. 4 ## analize the exponents. exx:= -rule[p]; for i from 1 by 1 to p-1 do exx:= exx +rule[p+i]- rule[i]; od: if type(exx,integer) then if exx > (p-2) then print(`Error: g[0],...,g[p-2] are not bases. You need a manual operation.`); RETURN(exx); fi: fi: if Verbose then print(`The exponent is `,exx); fi: L:=expand(x^(-1)*deltaToD(hg1(p-1) - x*hg2(p),t)); L:=expand(subs(getSubs(rule,p),L)); L1:=getL1(param,rule,p,+1); L1down:=getL1(param,rule,p,-1); L:=subs(x=x+1,L); L1:=subs(x=x+1,L1); L1down:=subs(x=x+1,L1down); degL1:=degree(L1[3],Dx); degL1down:=degree(L1down[3],Dx); tmp:=getBestDecomposition(r,degL1,degL1down); rup:=tmp[1]; rdown:=tmp[2]; if Verbose then print(`Our decompsition is `,tmp); fi: clist:=compose0(param,L1,rup); if invalid(clist) then print(`Error: H.G. equation is reducible. You need a manual operation.`); RETURN(0); fi; clist:=[op(myReverse(clist)),op(compose0down(param,L1down,rdown))]; if invalid(clist) then print(`Error: H.G. equation is reducible. You need a manual operation.`); RETURN(0); fi; subsrule:= getRule(sequence(L,max(rup*degL1,rdown*degL1down))); if Verbose then print(`We have gotten the rule.`); fi: mat:= getMat(base,r,clist,subsrule); if Verbose then print(`We have obtained the matrix to be solved.`); fi: mat:= factorout(mat); alist:=[]; for i from 1 by 1 to (r+1) do alist:=[op(alist),factor(clist[i][2]/((clist[i][1])*mat[2][i]))]; od: mat:=mat[1]; RETURN([getDifferenceOp(rup,rdown),alist,linsolve(mat,zeroVec(base+1))]); end: