share/contrib/Zeilberger/zeilberger_algorithm.mac

   1 /* Zeilberger's algorithm        */
   2
   3 /* RISC Institute, Linz, Austria */
   4 /* by Fabrizio Caruso            */
   5
   6
   7
   8 /* Routine that generates a parametrized function */
   9 parShiftAux(f,n,parName,degree,count) :=
  10   if count = degree then
  11      parName[count] * subst(n+count,n,f)
  12   else
  13      parName[count]*subst(n+count,n,f) +
  14      parShiftAux(f,n,parName,degree,count+1);
  15
  16
  17 /* Short name macro for parShift */
  18 parShift(f,n,parName,degree) ::=
  19   buildq([f,n,parName,degree],
  20          parShiftAux(f,n,parName,degree,0));
  21
  22
  23
  24 /* It normalizes a solution by clearing the denominators */
  25 normalizeGZ(Gsol,Zsol) :=
  26   block([i,lcm,Gres,Zres,lcm],
  27   lcm : 1,
  28   for i : 1 thru length(Gsol) do
  29     (
  30     lcm : lcm(lcm,denom(Gsol[i]))
  31     ),
  32   for i : 1 thru length(Zsol) do
  33     (
  34     lcm : lcm(lcm,denom(Zsol[i]))
  35
  36     ),
  37
  38   Gres : [],
  39   Zres : [],
  40
  41   for i : 1 thru length(Gsol) do
  42      Gres : endcons(factor(Gsol[i] * lcm), Gres),
  43   for i : 1 thru length(Zsol) do
  44      Zres : endcons(factor(Zsol[i] * lcm), Zres),
  45
  46   return([Gres,Zres])
  47   );
  48
  49
  50
  51 /* Constant Factor correction for the Zeilberger solutions */
  52 constFactCorr(sol,constFact,var) :=
  53   block([num,den,res,i,j,numCorr,denCorr],
  54      num : first(constFact),
  55      den : second(constFact),
  56      res : [],
  57
  58
  59 /* Remark: length(sol)-1 = ZAnsatzDegree */
  60
  61      for i : 1 thru length(sol) do
  62         (
  63         numCorr : 1,
  64         for j : 1 thru i-1 do
  65            numCorr : numCorr * subst(var+j-1,var,num),
  66
  67         denCorr : 1,
  68         for j : 1 thru length(sol)-i do
  69            denCorr : denCorr * subst(var+length(sol)-1-j,var,den),
  70
  71         res : endcons(factor(sol[i]/numCorr/denCorr),res)
  72
  73         ),
  74
  75      return(res)
  76   );
  77
  78
  79
  80 /* Separates the Gosper and Zeilberger Ansatz coefficients */
  81 separateGZVerboseOpt(Gres,GLength,ZLength,mode) :=
  82   block(
  83   [GList,ZList],
  84
  85   GList:[],
  86   ZList:[],
  87   if mode>= verbose then
  88     (
  89     print("Gres : " , Gres),
  90     print("GLength : " , GLength),
  91     print("ZLength : " , ZLength)
  92     ),
  93   for i : 1 thru GLength do
  94     (
  95     if mode>=verbose then
  96        print("Inserting ", Gres[i], " in the Gosper-list"),
  97     GList : endcons(rhs(Gres[i]),GList)
  98     ),
  99
 100   for i : GLength+1 thru GLength+ZLength do
 101     (
 102     if mode>=verbose then
 103        print("Inserting ", Gres[i], " in the Zeilberger-list"),
 104     ZList : endcons(rhs(Gres[i]),ZList)
 105     ),
 106
 107
 108   return([GList,ZList])
 109   );
 110
 111
 112
 113
 114 /* It creates the specific polynomial p */
 115 makeSpec(p) :=
 116    block(
 117    Spec : p,
 118    for i : 0 thru length(ZAnsatzCoeff)-1 do
 119      Spec : subst(factor(expand(ZAnsatzCoeff[i+1])),'_z[i],Spec),
 120    return(Spec)
 121    );
 122
 123 printParGosperSummary(p,sysSol,mode) :=
 124   block([],
 125    print("--------------------------"),
 126    print("Zeilberger-related summary"),
 127    print("(1) reconstructed p : ", p),
 128    print("(2) dimension of the solution : ", length(sysSol)),
 129    if mode>= verbose then
 130       print("(3) solution : ", sysSol)
 131    );
 132
 133
 134 /* Zeilberger's algorithm ("Parametrized" Gosper's algorithm) */
 135 parGosperVerboseOpt(f,k,n,ZAnsatzDegree,mode) :=
 136    block(
 137    [
 138    parF,           /* Parametrized function */
 139    parFRat,        /* Parametrized quotient function */
 140    nF,nFTrivial, numNF,denNF,constDenNF,
 141    GosperForm,     /* Gosper form */
 142    p,q,r,          /* Gosper form's polynomials */
 143    GAnsatzDegree,  /* Gosper Ansatz' degree */
 144    Gpoly,          /* Gosper polynomial */
 145    sysSol,         /* Solution to the Gosper system */
 146    normSysSol,     /* Normalized solution to the Gosper system */
 147    SpSysSol,       /* Specific solution to the Gosper system */
 148    GAnsatzCoeff,   /* Gosper Ansatz coefficients */
 149    ZAnsatzCoeff,   /* Zeilberger Ansatz coeffiecients */
 150    GZAnsatzCoeff,  /* Couple (GAnsatzCoeff,ZAnsatzCoeff) */
 151    GosperSol,      /* Solution to the Gosper equation */
 152    tlscope,        /* Solution to the telescoping problem */
 153    _g,             /* Gosper Ansatz coefficients */
 154    _z,              /* Zeilberger Ansatz coefficients */
 155    allSols,
 156    i,
 157    GosperSol,
 158    test_Gpoly,
 159    constantPart
 160    ],
 161
 162
 163
 164 if not(dependent(f,n)) and ZAnsatzDegree>0 and warnings then
 165   (
 166   print("---------------------------------------------"),
 167   print("WARNING!"),
 168   print("The input does not depend on ",n,"."),
 169   print("Maybe you should first try Gosper's algorithm."),
 170   print("---------------------------------------------")
 171   ),
 172
 173
 174 if ZAnsatzDegree = 0 then
 175   (
 176   if mode>=verbose then
 177     print("Invoking Gosper's algorithm"),
 178   GosperSol : GosperVerboseOpt(f,k,mode),
 179   if GosperSol=NO_HYP_SOL then
 180     return([])
 181   else
 182     if GosperSol = NON_HYPERGEOMETRIC then
 183       return([NON_PROPER_HYPERGEOMETRIC])
 184     else
 185       return([[GosperSol,[1]]])
 186   ),
 187
 188 nFTrivial : niceForm(f,n,'_z,ZAnsatzDegree,k),
 189
 190 if nFTrivial = [] then
 191   (
 192   if warnings then print(f, " is not hypergeometric in ", n),
 193   return([NON_PROPER_HYPERGEOMETRIC])
 194   ),
 195
 196 nF : nFTrivial[1],
 197
 198
 199 if mode >= verbose then
 200   (
 201   print("NICE FORM : ", nF),
 202   print("constant parts :  [", nFTrivial[2], ",", nFTrivial[3],"]")
 203   ),
 204
 205 numNF : num(nF),
 206 denNF : factor(denom(nF)),
 207
 208 if mode >= very_verbose then
 209   print("denNF : ", factor(denNF)),
 210
 211 if mode >= very_verbose then
 212    print("Considering now : ", factor(f/denNF)),
 213
 214 parFRat : shiftQuoOnlyHyp(factor(f/denNF),k), /* We temporarily forget about the factor numNF */
 215
 216 if parFRat = [] then
 217  (
 218  if warnings then print(f, " is not hypergeometric in ", k),
 219  return([NON_PROPER_HYPERGEOMETRIC])
 220  ),
 221
 222
 223 if mode >= verbose then
 224    print("--- Shift quotient computed! : " , parFRat),
 225
 226
 227 /* Computation of the Gosper form */
 228
 229    GosperForm: makeGosperFormVerboseOpt(parFRat,k,mode-1),
 230
 231
 232 gfP: takeP(GosperForm),
 233
 234 if mode >= very_verbose then
 235   print("gfP :", factor(gfP)),
 236
 237
 238 p: gfP * numNF, /* Now we remember that we "forgot" the factor numNF and compute the real p */
 239
 240 q: takeQ(GosperForm),
 241
 242 r: takeR(GosperForm),
 243
 244 if mode>= verbose then
 245   (
 246   print(" p : ", p),
 247   print(" q : ", q),
 248   print(" r : ", r)
 249   ),
 250
 251
 252 /* Solution of the Gosper equation */
 253 GAnsatzDegree: AnsatzDegreeVerboseOpt(p,q,r,k,mode-1),
 254
 255
 256   Gpoly : 0,
 257
 258   if integerp(GAnsatzDegree) and GAnsatzDegree >= 0 then
 259     Gpoly: GosperPoly(p,q,r,k,'_g,GAnsatzDegree),
 260
 261
 262 if mode >= very_verbose then
 263   print("Gosper's equation : ", Gpoly),
 264
 265
 266 if Gpoly = 0 then
 267   return([]),
 268
 269
 270
 271 /* --------------------------------------------------------------------- */
 272 /* Numerical test                                                        */
 273 if mod_test and (ZAnsatzDegree<=mod_threshold) then
 274   (
 275   modulus : mod_big_prime,
 276   test_Gpoly : polymod(subst(ev_point,n,Gpoly),mod_big_prime),
 277
 278
 279
 280   if test_Gpoly = 0 then
 281     return([]),
 282
 283   if mode>= very_verbose then
 284     print("modulus : ", modulus),
 285   sysSol : SolveZSysVerboseOpt(test_Gpoly,
 286                               '_g,GAnsatzDegree+1,'_z,ZAnsatzDegree+1,k,
 287                               modular_linear_solver,mode-1),
 288
 289   modulus : false,
 290   if mode>= very_verbose then
 291      print("numerical sysSol : ", sysSol),
 292   if sysSol = [] then
 293     (
 294     if mode >= verbose then
 295       print("no solution in the numerical test!"),
 296     return([])
 297     )
 298   else
 299     if mode >= verbose then
 300        (
 301        print("solution found in the numerical test!"),
 302        print("sysSol : ", sysSol)
 303        )
 304   ),
 305 /* End of numerical test                                                 */
 306 /* --------------------------------------------------------------------- */
 307
 308
 309 sysSol : SolveZSysVerboseOpt(Gpoly,'_g,GAnsatzDegree+1,'_z,ZAnsatzDegree+1,k,
 310                              linear_solver,mode-1),
 311
 312 if mode>= very_verbose then
 313   (
 314   print("sysSol: "),
 315   print(sysSol),
 316   print("%rnum_list : ", %rnum_list)
 317   ),
 318
 319 if warnings and length(sysSol)>1 then
 320   (
 321   print("------------------------------------------------------------"),
 322   print("WARNING!"),
 323   print("Multiple solutions have been found."),
 324   print("A lower order might probably work and yield fewer solutions."),
 325   print("------------------------------------------------------------")
 326   ),
 327
 328 SpSysSol : sysSol, /* non-sense? */
 329
 330 allSols : [],
 331 for i: 1 thru length(sysSol) do
 332   (
 333   if mode>=very_verbose then
 334     (
 335     print("Solution no. ", i),
 336     print(sysSol[i])
 337     ),
 338   SpSysSol : sysSol[i],
 339
 340
 341   GZAnsatzCoeff:separateGZVerboseOpt(SpSysSol,GAnsatzDegree+1,ZAnsatzDegree+1,
 342                                      mode-1),
 343
 344   GAnsatzCoeff: GZAnsatzCoeff[1],
 345   ZAnsatzCoeff: GZAnsatzCoeff[2],
 346
 347   if mode>=very_verbose then
 348     (
 349     print("GAnsatzCoeff : ", GAnsatzCoeff),
 350     print("ZAnsatzCoeff : ", ZAnsatzCoeff)
 351     ),
 352
 353
 354   ZAnsatzCoeff : constFactCorr(ZAnsatzCoeff,[nFTrivial[2],nFTrivial[3]],n),
 355
 356 if mode>=very_verbose then
 357   print("correction : ",[nFTrivial[2],nFTrivial[3]]),
 358
 359
 360
 361   constDenNF : product(subst(n+i-1,n,nFTrivial[3]),i,1,ZAnsatzDegree),
 362
 363
 364   if mode>=very_verbose then
 365     (
 366     print("constDenNF : ", constDenNF),
 367     print("CORRECTED ZAnsatzCoeff : ", ZAnsatzCoeff)
 368     ),
 369
 370   GZAnsatzCoeff : normalizeGZ(GAnsatzCoeff,ZAnsatzCoeff),
 371   GAnsatzCoeff : GZAnsatzCoeff[1],
 372
 373   ZAnsatzCoeff : GZAnsatzCoeff[2],
 374
 375
 376   if mode>=extra then
 377     (
 378     print(""),
 379     SolDeg : getDegrees(SpSysSol,GAnsatzDegree+1,ZAnsatzDegree+1,n),
 380     print("Degrees of the Gosper parameters"),
 381     prettyPrint(SolDeg[1]),
 382     print(""),
 383     print("Degrees of the Zeilberger parameters"),
 384     prettyPrint(SolDeg[2])
 385     ),
 386
 387
 388
 389 /*  GSol: factor(expand(sol2Poly(GAnsatzCoeff,k))), */
 390   GSol : sol2Poly(GAnsatzCoeff,k),
 391   if mode>= verbose then
 392     print("GSol : ", GSol),
 393   if Gsol = 0 then
 394      (
 395      grat:0,
 396      tlscope:0
 397      )
 398   else
 399      tlscope : factor(GSol * r / gfP / denNF / constDenNF),
 400
 401
 402
 403  if trivial_solutions or (not(tlscope=0) and
 404     not(apply("and",map(lambda([x],is(x=0)),ZAnsatzCoeff)))) then
 405
 406
 407    if simplified_output then
 408      allSols : endcons([ratsimp(tlscope),ZAnsatzCoeff],allSols)
 409    else
 410      allSols : endcons([tlscope,ZAnsatzCoeff],allSols)
 411  else
 412    if warnings then
 413      if tlscope = 0 then
 414        (
 415        print("----------------------------------------"),
 416        print("WARNING!"),
 417        print("The algorithm has discarted a trivial solution"),
 418        print("with zero certificate."),
 419
 420        if mode>= very_verbose then
 421         print("The coefficients of the recurrence are: ", ZAnsatzCoeff),
 422        print("----------------------------------------")
 423        )
 424      else
 425       (
 426       print("-------------------------------------------"),
 427       print("WARNING!"),
 428       print("The algorithm has discarted a trivial solution"),
 429       print("whose recurrence only has zero coefficients."),
 430
 431       if mode>= very_verbose then
 432          (
 433          print("The certificate is: "),
 434          print(tlscope)
 435          ),
 436       print("-------------------------------------------")
 437       )
 438  ),
 439
 440 if mode>=summary then
 441      (
 442      printGosperSummary(parFRat,GosperForm,GAnsatzDegree,Gpoly,GSol,mode-1),
 443      printParGosperSummary(p,sysSol,mode)
 444      ),
 445
 446   return(allSols)
 447
 448 );
 449
 450
 451
 452 /* It builds the l.h.s. of the recursion of the summmand */
 453 /* out of the solution of "parGosper" */
 454 makeRecursionAux(f,ZCoeff,n,count) :=
 455   if count = 1 then
 456      f * ZCoeff[1]
 457   else
 458      f * ZCoeff[1] +
 459      makeRecursionAux(subst(n+1,n,f),Rest(ZCoeff),n,count-1) ;
 460
 461
 462 /* Front-end macro of makeRecursionAux */
 463 makeRecursion(f,ZCoeff,n) ::=
 464   buildq([f,ZCoeff,n],
 465   makeRecursionAux(f,ZCoeff,n,length(ZCoeff))
 466         );
 467
 468
 469 UnEvMakeRecursion(f,n,ord) :=
 470   block(
 471   rec : concat(concat("a[",0,"]"),"f(n,k)"),
 472   for i : 1 thru ord do
 473      rec : concat(rec, "+", concat(concat("a[",i,"]"),concat("f(",string(n+i),",k)"))),
 474   rec
 475   );
 476
 477 /* It builds the l.h.s. of the recursion of the summmand */
 478 /* out of the solution of "parGosper" */
 479 makeRecursionOpAux(f,ZCoeff,n,count,upperBound) :=
 480   if count = upperBound then
 481      first(ZCoeff) * eta[n]^count ( f )
 482   else
 483      first(ZCoeff) * eta[n]^count ( f ) +
 484      makeRecursionOpAux(f,Rest(ZCoeff),n,count+1,upperBound) ;
 485
 486
 487 /* Front-end macro of makeRecursionAux */
 488 makeRecursionOp(f,ZCoeff,n) :=
 489    block(
 490    count : 1,
 491    makeRecursionOpAux(f,ZCoeff,n,count,length(ZCoeff))
 492    ) ;
 493
 494
 495 printZeilbergerSummary(f,parGRes) :=
 496   block([i],
 497
 498      print("---------------------------"),
 499      print("Iterated-Zeilberger summary"),
 500      print(UnEvMakeRecursion(f,n,length(parGRes[1][2])-1)," = ",
 501       Delta[k]("cert(n,k)"*"f(n,k)")),
 502
 503      print("where"),
 504      print("f(n,k) = ", f),
 505      print("and"),
 506      print("cert(n,k) = ", parGRes[1][1]),
 507      print("and"),
 508
 509
 510      for i : 0 thru length(parGRes[1][2])-1 do
 511       print("a[",i,"] = ", parGRes[1][2][i+1])
 512      );
 513
 514
 515 /* Front-end of the Zeilberger's algorithm with verbosity option */
 516 ZeilbergerVerboseOpt(f,k,n,mode) :=
 517    block(
 518    [
 519    start,
 520    ord,
 521    parGRes
 522    ],
 523
 524    if Gosper_in_Zeilberger then
 525       start : 0
 526    else
 527       (
 528       start : 1,
 529       if warnings and not(dependent(f,n)) then
 530         (
 531         print("-------------------------------------------------"),
 532         print("WARNING!                             "),
 533         print(f, " does not depend on ", n,"."),
 534         print("The flag Gosper_in_Zeilberger is set to false."),
 535         print("You might find a solution with Gosper's algorithm"),
 536         print("-------------------------------------------------")
 537         )
 538       ),
 539
 540
 541
 542    parGRes : [],
 543
 544    if mode>= verbose then
 545       print( "max_ord : " , max_ord ) ,
 546
 547    for ord : start step 1 while (parGRes = []) and ( ord <= max_ord ) do
 548       (
 549       if mode>= verbose then
 550         (
 551         print("----------------------------------------"),
 552         print("Searching a recurrence with order : " , ord )
 553         ),
 554
 555       parGRes : parGosperVerboseOpt(f,k,n,ord,mode-1),
 556       if mode >= verbose then
 557         (
 558         print("Gosper result with order : ", ord, " is : "),
 559         print(parGRes)
 560         )
 561       ),
 562
 563    if mode >= summary then
 564      printZeilbergerSummary(f,parGRes),
 565
 566    return(parGRes)
 567
 568    ) ;
 569
 570
 571 /* Non verbose version of Zeilberger's algorithm */
 572 ZeilbergerSummary([args]):=
 573  block([f,k,n],
 574  f:first(args),
 575  k:second(args),
 576  n:third(args),
 577  if length(args)=3 then
 578    ZeilbergerVerboseOpt(f,k,n,summary)
 579  else
 580    parGosperVerboseOpt(f,k,n,fourth(args),summary)
 581  );
 582
 583
 584 /* Non verbose version of Zeilberger's algorithm */
 585 ZeilbergerVerbose([args]):=
 586  block([f,k,n],
 587  f:first(args),
 588  k:second(args),
 589  n:third(args),
 590  if length(args)=3 then
 591    ZeilbergerVerboseOpt(f,k,n,verbose)
 592  else
 593    parGosperVerboseOpt(f,k,n,fourth(args),verbose)
 594  );
 595
 596
 597 /* Non verbose version of Zeilberger's algorithm */
 598 ZeilbergerVeryVerbose([args]):=
 599  block([f,k,n],
 600  f:first(args),
 601  k:second(args),
 602  n:third(args),
 603  if length(args)=3 then
 604    ZeilbergerVerboseOpt(f,k,n,very_verbose)
 605  else
 606    parGosperVerboseOpt(f,k,n,fourth(args),very_verbose)
 607  );
 608
 609
 610 /* Non verbose version of Zeilberger's algorithm */
 611 ZeilbergerExtra([args]):=
 612  block([f,k,n],
 613  f:first(args),
 614  k:second(args),
 615  n:third(args),
 616  if length(args)=3 then
 617    ZeilbergerVerboseOpt(f,k,n,extra)
 618  else
 619    parGosperVerboseOpt(f,k,n,fourth(args),extra)
 620  );
 621
 622
 623 /* Non verbose version of Zeilberger's algorithm */
 624 Zeilberger([args]):=
 625  block([f,k,n],
 626  f:first(args),
 627  k:second(args),
 628  n:third(args),
 629  if length(args)=3 then
 630    ZeilbergerVerboseOpt(f,k,n,nonverbose)
 631  else
 632    parGosperVerboseOpt(f,k,n,fourth(args),nonverbose)
 633  );
 634
 635
 636 /* Short name macro for verbose parGosper */
 637 parGosperVerbose(f,k,n,degree)::=
 638    buildq([f,k,n,degree],
 639           parGosperVerboseOpt(f,k,n,degree,verbose));
 640
 641
 642 /* Short name macro for very verbose parGosper */
 643 parGosperVeryVerbose(f,k,n,degree)::=
 644    buildq([f,k,n,degree],
 645           parGosperVerboseOpt(f,k,n,degree,very_verbose));
 646
 647
 648 /* Short name macro for the linear system version of parGosper */
 649 parGosperExtra(f,k,n,degree) ::=
 650    buildq([f,k,n,degree],
 651           parGosperVerboseOpt(f,k,n,degree,extra));
 652
 653 /* Short name macro for the summary version of parGosper */
 654 parGosperSummary(f,k,n,degree) ::=
 655    buildq([f,k,n,degree],
 656           parGosperVerboseOpt(f,k,n,degree,summary));
 657
 658
 659 /* Short name for the non-verbose version of parGosper */
 660 parGosperNonVerbose(f,k,n,degree) ::=
 661    buildq([f,k,n,degree],
 662           parGosperVerboseOpt(f,k,n,degree,nonverbose));
 663
 664 parGosper(f,k,n,order) ::=
 665    buildq([f,k,n,order],
 666           parGosperVerboseOpt(f,k,n,order,nonverbose));
 667
 668 /* Pretty print of a list */
 669 prettyPrint(list) :=
 670    block(
 671    len : length(list),
 672    for i : 1 step 1 thru len do
 673       list[i]
 674    );
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697