1
/* Gosper's algorithm */
3
/* RISC Institute, Linz, Austria */
4
/* by Fabrizio Caruso */
8
/* makeGosperForm2.macsyma */
9
/* poly2quint.macsyma */
10
/* GosperEq2.macsyma */
11
/* AlgUtil2.macsyma */
15
/* PREPROCESSING ROUTINES */
18
/* Preprocessing function "fRat" */
19
simpleFRatSlow(f,k) :=
21
factor((subst(k+1,k,f)/f)));
23
simpleFRatFast(f,k) :=
26
simpleFRat(f,k) := simpleFRatFast(f,k);
31
factor(expand(num(xthru(minfactorial(makefact(simpleFratSlow(f,k))))))/
32
expand(denom(xthru(minfactorial(makefact(simpleFratSlow(f,k))))))));
36
factor(expand(num(xthru(minfactorial(makefact(simpleFratFast(f,k))))))/
37
expand(denom(xthru(minfactorial(makefact(simpleFratFast(f,k)))))));
40
factor(expand(num(xthru(minfactorial(makefact(simpleFratFast(f,k))))))/
41
expand(denom(xthru(minfactorial(makefact(simpleFratFast(f,k)))))));
45
fRat(f,k) := fRatFast(f,k);
48
/* Postprocessing function that computes the rational function "certificate" */
49
telescope(Gsol,p,r,k) :=
52
/* Postprocessing function that computes the rational function "certificate" */
53
telescopeSimplified(Gsol,p,r,k) :=
57
/* Postprocessing function that computes the rational function "certificate" */
58
/* This is correct because we are telescoping w.r.t. a linear operator applied */
59
/* to the original function f , i.e. lin.op.(f) = (numNF/denNF)*f */
60
/* We telescoping solution is (r/p)*Gsol*(numNF/denNF)*f = */
61
/* = (r/gfP*numNF)*Gsol*(numNF/denNF) = (r*GSol)/(gfP*denNF) */
62
zTelescope(GSol,gfP,r,denNF) :=
65
/* Postprocessing function that computes the rational function "certificate" */
66
zTelescopeSimplified(GSol,gfP,r,denNF) :=
67
r*Gsol/(gfP*denNF); /* BOGUS */
70
printGosperSummary(fRat,GosperForm,AnsatzDegree,Gpoly,GosperSol,mode) :=
72
print("--------------------------------"),
73
print("Gosper-related summary"),
74
print("(1) fRat : " , fRat),
75
print("(2) Gosper form : ", GosperForm),
76
print("(3) AnsatzDegree : ", AnsatzDegree),
77
print("(4) Gosper polynomial : " , Gpoly),
78
if mode>= verbose then
79
print("(5) Gosper solution : ", GosperSol)
84
/* GOSPER'S ALGORITHM - Main routine - verbose */
85
GosperVerboseOpt(f,k,mode) :=
100
/* Computation of the shift quotient */
101
fRat : shiftQuoOnlyHyp(factor(f),k), /* factor might not be necessary */
104
print(f, " is not hypergeometric in ", k),
105
return(NON_HYPERGEOMETRIC)
108
if mode>=verbose then
109
print("Shift quotient computed! : " , fRat),
112
/* Computation of the Gosper form */
113
GosperForm: makeGosperFormVerboseOpt(factor(fRat),k,mode-1), /* factor ? */
115
if GosperForm = NON_HYPERGEOMETRIC then
117
print(f, " is not proper hypergeometric in ", k),
118
return(NON_HYPERGEOMETRIC)
121
p: takeP(GosperForm),
122
q: takeQ(GosperForm),
123
r: takeR(GosperForm),
124
if mode>=verbose then
125
print("p,q,r have been extracted! p: ", p, ", q: ", q, ", r: ", r),
127
/* Solution of the Gosper equation */
128
AnsatzDegree: AnsatzDegreeVerboseOpt(p,q,r,k,mode-1),
130
if mode>=verbose then
131
print("Ansatz degree computed ! : " , AnsatzDegree),
135
if AnsatzDegree >= 0 then
136
Gpoly: expand(GosperPoly(p,q,r,k,"c",AnsatzDegree)),
137
if mode>=verbose then
138
print("Gosper polynomial created! : " , Gpoly),
142
sysSol: SolveSysVerboseOpt(Gpoly,"c",AnsatzDegree+1,k,mode-1),
144
if mode>=very_verbose then
145
print("System of equations solved! : " , sysSol),
146
GosperSol: Gsol2Poly(first(sysSol),k),
147
if mode>=verbose then
148
print("Gosper solution created! : " , GosperSol),
150
if mode>=summary then
151
printGosperSummary(fRat,GosperForm,AnsatzDegree,Gpoly,GosperSol,mode-1),
154
/* Postprocessing : building the sol. to the telesc. probl. */
155
if GosperSol = 0 then
158
tlscope: telescope(GosperSol,p,r,k),
161
if simplified_output then
162
return(ratsimp(tlscope))
168
/* Macro for the normal verbose Gosper algorithm */
169
GosperVerbose(f,k) ::=
171
GosperVerboseOpt(f,k,verbose));
173
/* Macro for the very verbose Gosper algorithm */
174
GosperVeryVerbose(f,k) ::=
176
GosperVerboseOpt(f,k,very_verbose));
181
GosperVerboseOpt(f,k,extra));
183
/* Macro for the summary version of the Gosper algorithm */
184
GosperSummary(f,k) ::=
186
GosperVerboseOpt(f,k,summary));
189
/* Short name for the non verbose mode */
192
GosperVerboseOpt(f,k,nonverbose));
194
/* ----------------------------------------------------------------------- */
195
/* Anti-difference : It finds the anti-difference of a hypergeometric term */
197
AntiDifferenceVerboseOpt(f,k,mode) :=
199
gSol : GosperVerboseOpt(f,k,mode),
200
if gSol=NO_HYP_SOL then
201
return(NO_HYP_ANTIDIFFERENCE)
203
if gSol=NON_HYPERGEOMETRIC then
204
return(NON_HYPERGEOMETRIC)
206
if simplified_output then
207
return(factor(f*gSol))
213
AntiDifference(f,k) :=
214
AntiDifferenceVerboseOpt(f,k,nonverbose);
216
AntiDifferenceSummary(f,k) :=
217
AntiDifferenceVerboseOpt(f,k,summary);
219
AntiDifferenceVerbose(f,k) :=
220
AntiDifferenceVerboseOpt(f,k,verbose);
222
AntiDifferenceVeryVerbose(f,k) :=
223
AntiDifferenceVerboseOpt(f,k,very_verbose);
225
AntiDifferenceExtra(f,k) :=
226
AntiDifferenceVerboseOpt(f,k,extra);
229
/* ------------------------------------------------------------------ */
230
/* Gosper Sum : It sums Gosper-summable hypergeometric terms */
232
/* Gosper summation function - Summary mode */
233
GosperSumVerboseOpt(f,k,a,b,mode) :=
237
antiDiff:AntiDifferenceVerboseOpt(f,k,mode),
238
if antiDiff=NO_HYP_ANTIDIFFERENCE then
239
return(NON_GOSPER_SUMMABLE)
241
if antiDiff=NON_HYPERGEOMETRIC then
242
return(NON_HYPERGEOMETRIC)
245
if mode>= verbose then
246
print("Anti-difference : ", antiDiff),
247
return(subst(b+1,k,antiDiff)-subst(a,k,antiDiff))
251
/* Gosper summation function */
252
GosperSum(f,k,a,b) :=
253
GosperSumVerboseOpt(f,k,a,b,nonverbose);
256
/* Gosper summation function - Summary mode */
257
GosperSumSummary(f,k,a,b) :=
258
GosperSumVerboseOpt(f,k,a,b,summary);
260
GosperSumVerbose(f,k,a,b) :=
261
GosperSumVerboseOpt(f,k,a,b,verbose);
263
GosperSumVeryVerbose(f,k,a,b) :=
264
GosperSumVerboseOpt(f,k,a,b,very_verbose);
266
GosperSumExtra(f,k,a,b) :=
267
GosperSumVerboseOpt(f,k,a,b,extra);
added
removed
share/contrib/Zeilberger/Gosper.mac
