share/contrib/Zeilberger/Gosper.mac

   1 /* Gosper's algorithm        */
   2
   3 /* RISC Institute, Linz, Austria */
   4 /* by Fabrizio Caruso            */
   5
   6
   7 /* Dependencies:             */
   8 /* makeGosperForm2.macsyma   */
   9 /*      poly2quint.macsyma   */
  10 /* GosperEq2.macsyma         */
  11 /*      AlgUtil2.macsyma     */
  12
  13
  14
  15 /* PREPROCESSING ROUTINES */
  16
  17
  18 /* Preprocessing function "fRat" */
  19 simpleFRatSlow(f,k) :=
  20   minfactorial(
  21   factor((subst(k+1,k,f)/f)));
  22
  23 simpleFRatFast(f,k) :=
  24   subst(k+1,k,f)/f;
  25
  26 simpleFRat(f,k) := simpleFRatFast(f,k);
  27
  28
  29 fRatSlow(f,k) :=
  30   minfactorial(
  31   factor(expand(num(xthru(minfactorial(makefact(simpleFratSlow(f,k))))))/
  32          expand(denom(xthru(minfactorial(makefact(simpleFratSlow(f,k))))))));
  33
  34
  35 fratFast(f,k) :=
  36   factor(expand(num(xthru(minfactorial(makefact(simpleFratFast(f,k))))))/
  37          expand(denom(xthru(minfactorial(makefact(simpleFratFast(f,k)))))));
  38
  39 fratVeryFast(f,k) :=
  40   factor(expand(num(xthru(minfactorial(makefact(simpleFratFast(f,k))))))/
  41          expand(denom(xthru(minfactorial(makefact(simpleFratFast(f,k)))))));
  42
  43
  44
  45 fRat(f,k) := fRatFast(f,k);
  46
  47
  48 /* Postprocessing function that computes the rational function "certificate" */
  49 telescope(Gsol,p,r,k) :=
  50    r*Gsol/p;
  51
  52 /* Postprocessing function that computes the rational function "certificate" */
  53 telescopeSimplified(Gsol,p,r,k) :=
  54    r*Gsol/factor(p);
  55
  56
  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) :=
  63    r*GSol/(gfP*denNF);
  64
  65 /* Postprocessing function that computes the rational function "certificate" */
  66 zTelescopeSimplified(GSol,gfP,r,denNF) :=
  67    r*Gsol/(gfP*denNF); /* BOGUS */
  68
  69
  70 printGosperSummary(fRat,GosperForm,AnsatzDegree,Gpoly,GosperSol,mode) :=
  71   block([],
  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)
  80      );
  81
  82
  83
  84 /* GOSPER'S ALGORITHM - Main routine - verbose */
  85 GosperVerboseOpt(f,k,mode) :=
  86   block(
  87   [fRat,
  88    signFlag,
  89    GosperForm,
  90    p,q,r,
  91    AnsatzDegree,
  92    Gpoly,
  93    sysSol,
  94    GosperSol,
  95    tlscope,
  96    c,
  97    constantPart],
  98
  99
 100
 101   /* Computation of the shift quotient */
 102   fRat : shiftQuoOnlyHyp(factor(f),k), /* factor might not be necessary */
 103   if fRat = [] then
 104     (
 105       if warnings then
 106         print(f, " is not hypergeometric in ", k),
 107     return(NON_HYPERGEOMETRIC)
 108     ),
 109
 110   if mode>=verbose then
 111      print("Shift quotient computed! : " , fRat),
 112
 113
 114   /* Computation of the Gosper form */
 115   GosperForm: makeGosperFormVerboseOpt(factor(fRat),k,mode-1), /* factor ? */
 116
 117   if GosperForm = NON_HYPERGEOMETRIC then
 118     (
 119     print(f, " is not proper hypergeometric in ", k),
 120     return(NON_HYPERGEOMETRIC)
 121     ),
 122
 123   p: takeP(GosperForm),
 124   q: takeQ(GosperForm),
 125   r: takeR(GosperForm),
 126   if mode>=verbose then
 127      print("p,q,r have been extracted! p: ", p, ", q: ", q, ", r: ", r),
 128
 129   /* Solution of the Gosper equation */
 130   AnsatzDegree: AnsatzDegreeVerboseOpt(p,q,r,k,mode-1),
 131
 132   if mode>=verbose then
 133      print("Ansatz degree computed ! : " , AnsatzDegree),
 134
 135   Gpoly : 0,
 136
 137   if AnsatzDegree >= 0 then
 138     Gpoly: expand(GosperPoly(p,q,r,k,c,AnsatzDegree)),
 139   if mode>=verbose then
 140      print("Gosper polynomial created! : " , Gpoly),
 141   if Gpoly = 0 then
 142      return(NO_HYP_SOL),
 143
 144   sysSol: SolveSysVerboseOpt(Gpoly,c,AnsatzDegree+1,k,mode-1),
 145
 146   if mode>=very_verbose then
 147      print("System of equations solved! :  " , sysSol),
 148   GosperSol: Gsol2Poly(first(sysSol),k),
 149   if mode>=verbose then
 150      print("Gosper solution created! : " , GosperSol),
 151
 152   if mode>=summary then
 153      printGosperSummary(fRat,GosperForm,AnsatzDegree,Gpoly,GosperSol,mode-1),
 154
 155
 156   /* Postprocessing : building the sol. to the telesc. probl. */
 157   if GosperSol = 0 then
 158     tlscope: NO_HYP_SOL
 159   else
 160     tlscope: telescope(GosperSol,p,r,k),
 161
 162
 163  if simplified_output then
 164     return(ratsimp(tlscope))
 165  else
 166     return(tlscope)
 167   );
 168
 169
 170 /* Macro for the normal verbose Gosper algorithm */
 171 GosperVerbose(f,k) ::=
 172   buildq([f,k],
 173          GosperVerboseOpt(f,k,verbose));
 174
 175 /* Macro for the very verbose Gosper algorithm */
 176 GosperVeryVerbose(f,k) ::=
 177   buildq([f,k],
 178          GosperVerboseOpt(f,k,very_verbose));
 179
 180
 181 GosperExtra(f,k) ::=
 182   buildq([f,k],
 183          GosperVerboseOpt(f,k,extra));
 184
 185 /* Macro for the summary version of the Gosper algorithm */
 186 GosperSummary(f,k) ::=
 187   buildq([f,k],
 188          GosperVerboseOpt(f,k,summary));
 189
 190
 191 /* Short name for the non verbose mode */
 192 Gosper(f,k) ::=
 193    buildq([f,k],
 194           GosperVerboseOpt(f,k,nonverbose));
 195
 196 /* ----------------------------------------------------------------------- */
 197 /* Anti-difference : It finds the anti-difference of a hypergeometric term */
 198
 199 AntiDifferenceVerboseOpt(f,k,mode) :=
 200   block([gSol],
 201   gSol : GosperVerboseOpt(f,k,mode),
 202   if gSol=NO_HYP_SOL then
 203     return(NO_HYP_ANTIDIFFERENCE)
 204   else
 205     if gSol=NON_HYPERGEOMETRIC then
 206       return(NON_HYPERGEOMETRIC)
 207     else
 208       if simplified_output then
 209          return(factor(f*gSol))
 210       else
 211          return(f*gSol)
 212   );
 213
 214
 215 AntiDifference(f,k) :=
 216   AntiDifferenceVerboseOpt(f,k,nonverbose);
 217
 218 AntiDifferenceSummary(f,k) :=
 219   AntiDifferenceVerboseOpt(f,k,summary);
 220
 221 AntiDifferenceVerbose(f,k) :=
 222   AntiDifferenceVerboseOpt(f,k,verbose);
 223
 224 AntiDifferenceVeryVerbose(f,k) :=
 225   AntiDifferenceVerboseOpt(f,k,very_verbose);
 226
 227 AntiDifferenceExtra(f,k) :=
 228   AntiDifferenceVerboseOpt(f,k,extra);
 229
 230
 231 /* ------------------------------------------------------------------ */
 232 /* Gosper Sum : It sums Gosper-summable hypergeometric terms          */
 233
 234 /* Gosper summation function - Summary mode */
 235 GosperSumVerboseOpt(f,k,a,b,mode) :=
 236   block(
 237         [antiDiff],
 238
 239         antiDiff:AntiDifferenceVerboseOpt(f,k,mode),
 240         if antiDiff=NO_HYP_ANTIDIFFERENCE then
 241            return(NON_GOSPER_SUMMABLE)
 242         else
 243           if antiDiff=NON_HYPERGEOMETRIC then
 244              return(NON_HYPERGEOMETRIC)
 245           else
 246             (
 247             if mode>= verbose then
 248               print("Anti-difference : ", antiDiff),
 249             return(subst(b+1,k,antiDiff)-subst(a,k,antiDiff))
 250             )
 251        );
 252
 253 /* Gosper summation function */
 254 GosperSum(f,k,a,b) :=
 255   GosperSumVerboseOpt(f,k,a,b,nonverbose);
 256
 257
 258 /* Gosper summation function - Summary mode */
 259 GosperSumSummary(f,k,a,b) :=
 260   GosperSumVerboseOpt(f,k,a,b,summary);
 261
 262 GosperSumVerbose(f,k,a,b) :=
 263   GosperSumVerboseOpt(f,k,a,b,verbose);
 264
 265 GosperSumVeryVerbose(f,k,a,b) :=
 266   GosperSumVerboseOpt(f,k,a,b,very_verbose);
 267
 268 GosperSumExtra(f,k,a,b) :=
 269   GosperSumVerboseOpt(f,k,a,b,extra);
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288