BeginPackage["NewZeil`"] Print["This is the Mathematica Package NewZeil"]; Print["version 0.99 , 12 April 2005"]; Print["written by Andrew Sills"]; Print["to accompany the paper"]; Print["`Sharp Upper Bounds for the Orders of the Recurrences Output by the Zeilberger and q-Zeilberger Algorithms'"]; Print["by Mohamud Mohammed and Doron Zeilberger"]; Print["J. Symbolic Comput., 39 (2005), 201-207. "]; Print[" "]; Print["Latest version available at http://www.math.rutgers.edu/~asills."]; Print["Please report bugs and/or send comments to "]; Print["asills at math dot rutgers dot edu"]; Print["For a listing of available functions, enter 'NewZeil[]'."]; NewZeil::usage = "NewZeil[] prints the list of functions in this package which are available to the user. The notation used in coding and documenting this package intentionally follows that of Mohammed and Zeilberger's paper `Sharp upper bounds...', JSC 39, 201-207 as closely as possible. " rf::usage = "rf[a,n] is an abbreviation for the standard Mathematica symbol Pochhammer[a,n]." Zeil::usage= "Zeil[H, z , POL, k, n, params], where H is the pure hypergeometric part of 'ProperHypergeometric' less the argument z (i.e. H(n,k)/z^k in the paper), z is the argument, POL is the polynomial part, k is the summation variable, n is the discrete parameter, params is the (optional) list of additional (possibly continuous) parameters; employs the simplified Zeilberger algorithm to { {Rec}, {Cert} }, where 'Rec' is a recurrence satisfied by Sum[ H z^k POL, {k,-Infinity,Infinity}] and 'Cert' is the rational function certificate R(n,k). " ; Begin["`Private`"] Verbose = False; NewZeil[]:=Module[{}, Print["User accessible functions available in NewZeil include"]; Print["rf[a,n], Zeil[H, z, POL, k, n, params ]."]; Print[" 'Verbose' is a global variable."]; Print[" Enter ? for help on individual functions."]; ]; rf[a_,n_]:=Module[{}, Pochhammer[a,n] ]; FindDeg[f_,g_,h_,k_]:= Module[{ degX, lam, A, B, d}, degX = Exponent[h,k] - Max[ Exponent[f,k], Exponent[g,k] ]; If[ Exponent[f,k] == Exponent[g,k] && Coefficient[f,k,Exponent[f,k]] == Coefficient[g,k,Exponent[g,k]], lam = Coefficient[ g, k, Exponent[g,k] ]; B = Coefficient[ g, k, Exponent[g,k]-1 ]; A = Coefficient[ f, k, Exponent[f,k]-1 ]; d = Factor[ Together[ (B-A)/lam ]]; If [ IntegerQ[d], degX = Max[ 1 + degX, d ], ++degX ] ]; Return[degX] ]; ZeilOrderL[H_, z_, POL_, k_, n_, L_Integer, params_List:{}]:= Module[{A,e,f,g,h,H1,Hbar,HbarRat,Hbarden,HbarOverH,Hnum,Hden,M,R,x,X, vars,eqns,soln,CleanUpSubs,CleanR,Rec,a,app,ap,j,m,s,b,c, trialvarsubs,trialsoln,trialeqns,commondenom,freevar}, e = "e"; x = "x"; H1 = ReplaceAll[ H, {Binomial[j_,m_] -> j!/m!/(j-m)!, Pochhammer[a_,j_] -> (a+j-1)!/(a-1)!} ]; If[ Verbose, Print["============================================"]]; If[ Verbose , Print["Attempting to find a recurrence of order ",L]]; H2 = ReplaceAll [H1, {Binomial[j_,m_] -> Pochhammer[1,j]/Pochhammer[1,m]/Pochhammer[1,j-m], Factorial[m_]-> Pochhammer[1,m]} ]; If[Verbose, Print["Here is the form of the pure hypergeometric portion of the summand being used: ",H1 z^k] ]; If[Verbose, Print["Here is the alternate form: ",H2 z^k]] ; Hnum = Numerator[H1]; Hden = Denominator[H1]; HbarDen = Hden/.n->n+L; Hbar = Hnum/HbarDen; (* h =FullSimplify[Factor[Expand[ Sum[Subscript[e,i] FullSimplify[((POL * H1) /.n->n+i) /Hbar], {i,0,L}] ]]]; *) h= Sum[ Subscript[e,i] * (POL/. n -> n+i) * ( (Numerator[H2]/.Pochhammer[app_.+b_. , ap_. n + a_. ]-> Pochhammer[b+app+ap n + a, i ap]) * (Denominator[H2]/. Pochhammer[app_.+b_. , ap_. n + a_.]-> Pochhammer[b+app + ap n + a, (L-i) ap]) ) /. Pochhammer[a_,b_]->1, {i,0,L}]; HbarRat = FullSimplify[(Hbar/.k->k+1)/Hbar]; HbarOverH = 1/((Denominator[H2]/. Pochhammer[app_. + b_., ap_. n + a_.]-> Pochhammer[1+ b + app + ap n + a -1, ap L]) /. Pochhammer[a_,b_] -> 1); f = z Numerator[HbarRat]; g = Denominator[HbarRat]; M = FindDeg[f,g,h,k]; X = Sum[ Subscript[x,i] k^i, {i,0,M}]; GosperEq = f (X/.k->k+1) - (g/.k->k-1) X - h; If[Verbose, Print["Hbar = ", Hbar]; Print["HbarRat=", HbarRat]; Print["HbarOverH=",HbarOverH]; Print["f = ", f]; Print["g = ", g]; Print["h = ", h]; Print["M = ", M]; Print["X = ", X]; Print["GosperLHS =", Collect[GosperEq,k] ] ]; vars = Join[ Table[ Subscript[e,i],{i,0,L}], Table[ Subscript[x,i],{i,0,M}] ]; eqns = Table[ Coefficient[GosperEq,k,i]==0, {i,0, Exponent[GosperEq,k]}]; eqns = Complement[eqns, {True}]; If[Verbose, Print["Homogeneous system to solve:", eqns] ]; If[Verbose, Print[ Table[ Coefficient[eqns[[i]][[1]], vars[[j]] ] , {i,Length[eqns]},{j,Length[vars]}]//MatrixForm, vars//MatrixForm ] ]; trialvarsubs = Table[ params[[i]]-> Random[Integer,{1,10}],{i,Length[params]} ]; trialvarsubs = Union[ {n->Random[Integer,{1,10}]}, trialvarsubs]; trialeqns = eqns/.trialvarsubs; If[z[[0]]==Symbol, trialeqns = trialeqns/.{z->Random[Integer,{1,10}]}]; trialsoln = Solve[trialeqns,vars]; If[Verbose, Print["Trial numeric solution", trialsoln]]; If[ Flatten[Union[ Table[ Subscript[e,i]/.trialsoln,{i,0,L} ]]]==={0}, Return[{{0},{0}}], Print["It is probable that solution of order ",L," exists."]; Print["Please be patient while I solve the full symbolic system."]; freevar=vars[[-1]]; j = -1; While[ ((freevar/.trialsoln)== {0} ) && ( -j < Length[vars] ), freevar = vars[[j]]; j--; ]; If[ -j >= Length[vars], Print["Something seems to have gone wrong."]; Print["Please send bug report to asills@math.rutgers.edu"]; Return[Null] ]; If[Verbose, Print["Free variable is ",freevar,"; will be set to 1."] ]; eqns = Union[ eqns, {freevar==1 }]; soln = FullSimplify[Solve[eqns,vars]]; If[Verbose, Print["Solutions to system:", soln] ]; commondenom = Apply[ PolynomialLCM, Denominator[vars/.soln] ]; commondenom = Apply[PolynomialLCM,Denominator[Flatten[vars/.soln]]]; If[Verbose, Print["Common denominator: ", commondenom] ]; R = Factor[( (g/.k->k-1) X HbarOverH/ POL )/. soln ]; If[Verbose, Print["R =", R] ]; Rec = (Sum[ commondenom Subscript[e,i] "F"[n+i,k], {i,0,L}]/.soln) == "G"[n,k+1]-"G"[n,k]; Return[Simplify[ { Rec, commondenom R}]] ] ]; Zeil[ F_, z_, POL_, k_, n_, params_List:{}]:= Module[{i,ret}, L=0; ret = { {0}, {0} }; While[ ret == { {0},{0} }, ret=ZeilOrderL[F,z,POL,k,n,L,params]; L++ ]; Return[ret] ]; End[] EndPackage[]