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],
  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    test_Gpoly,
 158    constantPart,
 159    GSol,
 160    gfP
 161    ],
 162
 163
 164
 165 if not(dependent(f,n)) and ZAnsatzDegree>0 and warnings then
 166   (
 167   print("---------------------------------------------"),
 168   print("WARNING!"),
 169   print("The input does not depend on ",n,"."),
 170   print("Maybe you should first try Gosper's algorithm."),
 171   print("---------------------------------------------")
 172   ),
 173
 174
 175 if ZAnsatzDegree = 0 then
 176   (
 177   if mode>=verbose then
 178     print("Invoking Gosper's algorithm"),
 179   GosperSol : GosperVerboseOpt(f,k,mode),
 180   if GosperSol=NO_HYP_SOL then
 181     return([])
 182   else
 183     if GosperSol = NON_HYPERGEOMETRIC then
 184       return([NON_PROPER_HYPERGEOMETRIC])
 185     else
 186       return([[GosperSol,[1]]])
 187   ),
 188
 189 nFTrivial : niceForm(f,n,'_z,ZAnsatzDegree,k),
 190
 191 if nFTrivial = [] then
 192   (
 193   if warnings then print(f, " is not hypergeometric in ", n),
 194   return([NON_PROPER_HYPERGEOMETRIC])
 195   ),
 196
 197 nF : nFTrivial[1],
 198
 199
 200 if mode >= verbose then
 201   (
 202   print("NICE FORM : ", nF),
 203   print("constant parts :  [", nFTrivial[2], ",", nFTrivial[3],"]")
 204   ),
 205
 206 numNF : num(nF),
 207 denNF : factor(denom(nF)),
 208
 209 if mode >= very_verbose then
 210   print("denNF : ", factor(denNF)),
 211
 212 if mode >= very_verbose then
 213    print("Considering now : ", factor(f/denNF)),
 214
 215 parFRat : shiftQuoOnlyHyp(factor(f/denNF),k), /* We temporarily forget about the factor numNF */
 216
 217 if parFRat = [] then
 218  (
 219  if warnings then print(f, " is not hypergeometric in ", k),
 220  return([NON_PROPER_HYPERGEOMETRIC])
 221  ),
 222
 223
 224 if mode >= verbose then
 225    print("--- Shift quotient computed! : " , parFRat),
 226
 227
 228 /* Computation of the Gosper form */
 229
 230    GosperForm: makeGosperFormVerboseOpt(parFRat,k,mode-1),
 231
 232
 233 gfP: takeP(GosperForm),
 234
 235 if mode >= very_verbose then
 236   print("gfP :", factor(gfP)),
 237
 238
 239 p: gfP * numNF, /* Now we remember that we "forgot" the factor numNF and compute the real p */
 240
 241 q: takeQ(GosperForm),
 242
 243 r: takeR(GosperForm),
 244
 245 if mode>= verbose then
 246   (
 247   print(" p : ", p),
 248   print(" q : ", q),
 249   print(" r : ", r)
 250   ),
 251
 252
 253 /* Solution of the Gosper equation */
 254 GAnsatzDegree: AnsatzDegreeVerboseOpt(p,q,r,k,mode-1),
 255
 256
 257   Gpoly : 0,
 258
 259   if integerp(GAnsatzDegree) and GAnsatzDegree >= 0 then
 260     Gpoly: GosperPoly(p,q,r,k,'_g,GAnsatzDegree),
 261
 262
 263 if mode >= very_verbose then
 264   print("Gosper's equation : ", Gpoly),
 265
 266
 267 if Gpoly = 0 then
 268   return([]),
 269
 270
 271
 272 /* --------------------------------------------------------------------- */
 273 /* Numerical test                                                        */
 274 if mod_test and (ZAnsatzDegree<=mod_threshold) then
 275   (
 276   modulus : mod_big_prime,
 277   test_Gpoly : polymod(subst(ev_point,n,Gpoly),mod_big_prime),
 278
 279
 280
 281   if test_Gpoly = 0 then
 282     return([]),
 283
 284   if mode>= very_verbose then
 285     print("modulus : ", modulus),
 286   sysSol : SolveZSysVerboseOpt(test_Gpoly,
 287                               '_g,GAnsatzDegree+1,'_z,ZAnsatzDegree+1,k,
 288                               modular_linear_solver,mode-1),
 289
 290   modulus : false,
 291   if mode>= very_verbose then
 292      print("numerical sysSol : ", sysSol),
 293   if sysSol = [] then
 294     (
 295     if mode >= verbose then
 296       print("no solution in the numerical test!"),
 297     return([])
 298     )
 299   else
 300     if mode >= verbose then
 301        (
 302        print("solution found in the numerical test!"),
 303        print("sysSol : ", sysSol)
 304        )
 305   ),
 306 /* End of numerical test                                                 */
 307 /* --------------------------------------------------------------------- */
 308
 309
 310 sysSol : SolveZSysVerboseOpt(Gpoly,'_g,GAnsatzDegree+1,'_z,ZAnsatzDegree+1,k,
 311                              linear_solver,mode-1),
 312
 313 if mode>= very_verbose then
 314   (
 315   print("sysSol: "),
 316   print(sysSol),
 317   print("%rnum_list : ", %rnum_list)
 318   ),
 319
 320 if warnings and length(sysSol)>1 then
 321   (
 322   print("------------------------------------------------------------"),
 323   print("WARNING!"),
 324   print("Multiple solutions have been found."),
 325   print("A lower order might probably work and yield fewer solutions."),
 326   print("------------------------------------------------------------")
 327   ),
 328
 329 SpSysSol : sysSol, /* non-sense? */
 330
 331 allSols : [],
 332 for i: 1 thru length(sysSol) do
 333   (
 334   if mode>=very_verbose then
 335     (
 336     print("Solution no. ", i),
 337     print(sysSol[i])
 338     ),
 339   SpSysSol : sysSol[i],
 340
 341
 342   GZAnsatzCoeff:separateGZVerboseOpt(SpSysSol,GAnsatzDegree+1,ZAnsatzDegree+1,
 343                                      mode-1),
 344
 345   GAnsatzCoeff: GZAnsatzCoeff[1],
 346   ZAnsatzCoeff: GZAnsatzCoeff[2],
 347
 348   if mode>=very_verbose then
 349     (
 350     print("GAnsatzCoeff : ", GAnsatzCoeff),
 351     print("ZAnsatzCoeff : ", ZAnsatzCoeff)
 352     ),
 353
 354
 355   ZAnsatzCoeff : constFactCorr(ZAnsatzCoeff,[nFTrivial[2],nFTrivial[3]],n),
 356
 357 if mode>=very_verbose then
 358   print("correction : ",[nFTrivial[2],nFTrivial[3]]),
 359
 360
 361
 362   constDenNF : product(subst(n+i-1,n,nFTrivial[3]),i,1,ZAnsatzDegree),
 363
 364
 365   if mode>=very_verbose then
 366     (
 367     print("constDenNF : ", constDenNF),
 368     print("CORRECTED ZAnsatzCoeff : ", ZAnsatzCoeff)
 369     ),
 370
 371   GZAnsatzCoeff : normalizeGZ(GAnsatzCoeff,ZAnsatzCoeff),
 372   GAnsatzCoeff : GZAnsatzCoeff[1],
 373
 374   ZAnsatzCoeff : GZAnsatzCoeff[2],
 375
 376
 377   if mode>=extra then
 378     (
 379     print(""),
 380     SolDeg : getDegrees(SpSysSol,GAnsatzDegree+1,ZAnsatzDegree+1,n),
 381     print("Degrees of the Gosper parameters"),
 382     prettyPrint(SolDeg[1]),
 383     print(""),
 384     print("Degrees of the Zeilberger parameters"),
 385     prettyPrint(SolDeg[2])
 386     ),
 387
 388
 389
 390 /*  GSol: factor(expand(sol2Poly(GAnsatzCoeff,k))), */
 391   GSol : sol2Poly(GAnsatzCoeff,k),
 392   if mode>= verbose then
 393     print("GSol : ", GSol),
 394   if Gsol = 0 then
 395      (
 396      grat:0,
 397      tlscope:0
 398      )
 399   else
 400      tlscope : factor(GSol * r / gfP / denNF / constDenNF),
 401
 402
 403
 404  if trivial_solutions or (not(tlscope=0) and
 405     not(apply("and",map(lambda([x],is(x=0)),ZAnsatzCoeff)))) then
 406
 407
 408    if simplified_output then
 409      allSols : endcons([ratsimp(tlscope),ZAnsatzCoeff],allSols)
 410    else
 411      allSols : endcons([tlscope,ZAnsatzCoeff],allSols)
 412  else
 413    if warnings then
 414      if tlscope = 0 then
 415        (
 416        print("----------------------------------------"),
 417        print("WARNING!"),
 418        print("The algorithm has discarted a trivial solution"),
 419        print("with zero certificate."),
 420
 421        if mode>= very_verbose then
 422         print("The coefficients of the recurrence are: ", ZAnsatzCoeff),
 423        print("----------------------------------------")
 424        )
 425      else
 426       (
 427       print("-------------------------------------------"),
 428       print("WARNING!"),
 429       print("The algorithm has discarted a trivial solution"),
 430       print("whose recurrence only has zero coefficients."),
 431
 432       if mode>= very_verbose then
 433          (
 434          print("The certificate is: "),
 435          print(tlscope)
 436          ),
 437       print("-------------------------------------------")
 438       )
 439  ),
 440
 441 if mode>=summary then
 442      (
 443      printGosperSummary(parFRat,GosperForm,GAnsatzDegree,Gpoly,GSol,mode-1),
 444      printParGosperSummary(p,sysSol,mode)
 445      ),
 446
 447   return(allSols)
 448
 449 );
 450
 451
 452
 453 /* It builds the l.h.s. of the recursion of the summmand */
 454 /* out of the solution of "parGosper" */
 455 makeRecursionAux(f,ZCoeff,n,count) :=
 456   if count = 1 then
 457      f * ZCoeff[1]
 458   else
 459      f * ZCoeff[1] +
 460      makeRecursionAux(subst(n+1,n,f),Rest(ZCoeff),n,count-1) ;
 461
 462
 463 /* Front-end macro of makeRecursionAux */
 464 makeRecursion(f,ZCoeff,n) ::=
 465   buildq([f,ZCoeff,n],
 466   makeRecursionAux(f,ZCoeff,n,length(ZCoeff))
 467         );
 468
 469
 470 UnEvMakeRecursion(f,n,ord) :=
 471   block(
 472   rec : concat(concat("a[",0,"]"),"f(n,k)"),
 473   for i : 1 thru ord do
 474      rec : concat(rec, "+", concat(concat("a[",i,"]"),concat("f(",string(n+i),",k)"))),
 475   rec
 476   );
 477
 478 /* It builds the l.h.s. of the recursion of the summmand */
 479 /* out of the solution of "parGosper" */
 480 makeRecursionOpAux(f,ZCoeff,n,count,upperBound) :=
 481   if count = upperBound then
 482      first(ZCoeff) * eta[n]^count ( f )
 483   else
 484      first(ZCoeff) * eta[n]^count ( f ) +
 485      makeRecursionOpAux(f,Rest(ZCoeff),n,count+1,upperBound) ;
 486
 487
 488 /* Front-end macro of makeRecursionAux */
 489 makeRecursionOp(f,ZCoeff,n) :=
 490    block(
 491    count : 1,
 492    makeRecursionOpAux(f,ZCoeff,n,count,length(ZCoeff))
 493    ) ;
 494
 495
 496 printZeilbergerSummary(f,parGRes) :=
 497   block([i],
 498
 499      print("---------------------------"),
 500      print("Iterated-Zeilberger summary"),
 501      print(UnEvMakeRecursion(f,n,length(parGRes[1][2])-1)," = ",
 502       Delta[k]("cert(n,k)"*"f(n,k)")),
 503
 504      print("where"),
 505      print("f(n,k) = ", f),
 506      print("and"),
 507      print("cert(n,k) = ", parGRes[1][1]),
 508      print("and"),
 509
 510
 511      for i : 0 thru length(parGRes[1][2])-1 do
 512       print("a[",i,"] = ", parGRes[1][2][i+1])
 513      );
 514
 515
 516 /* Front-end of the Zeilberger's algorithm with verbosity option */
 517 ZeilbergerVerboseOpt(f,k,n,mode) :=
 518    block(
 519    [
 520    start,
 521    ord,
 522    parGRes
 523    ],
 524
 525    if Gosper_in_Zeilberger then
 526       start : 0
 527    else
 528       (
 529       start : 1,
 530       if warnings and not(dependent(f,n)) then
 531         (
 532         print("-------------------------------------------------"),
 533         print("WARNING!                             "),
 534         print(f, " does not depend on ", n,"."),
 535         print("The flag Gosper_in_Zeilberger is set to false."),
 536         print("You might find a solution with Gosper's algorithm"),
 537         print("-------------------------------------------------")
 538         )
 539       ),
 540
 541
 542
 543    parGRes : [],
 544
 545    if mode>= verbose then
 546       print( "max_ord : " , max_ord ) ,
 547
 548    for ord : start step 1 while (parGRes = []) and ( ord <= max_ord ) do
 549       (
 550       if mode>= verbose then
 551         (
 552         print("----------------------------------------"),
 553         print("Searching a recurrence with order : " , ord )
 554         ),
 555
 556       parGRes : parGosperVerboseOpt(f,k,n,ord,mode-1),
 557       if mode >= verbose then
 558         (
 559         print("Gosper result with order : ", ord, " is : "),
 560         print(parGRes)
 561         )
 562       ),
 563
 564    if mode >= summary then
 565      printZeilbergerSummary(f,parGRes),
 566
 567    return(parGRes)
 568
 569    ) ;
 570
 571
 572 /* Non verbose version of Zeilberger's algorithm */
 573 ZeilbergerSummary([args]):=
 574  block([f,k,n],
 575  f:first(args),
 576  k:second(args),
 577  n:third(args),
 578  if length(args)=3 then
 579    ZeilbergerVerboseOpt(f,k,n,summary)
 580  else
 581    parGosperVerboseOpt(f,k,n,fourth(args),summary)
 582  );
 583
 584
 585 /* Non verbose version of Zeilberger's algorithm */
 586 ZeilbergerVerbose([args]):=
 587  block([f,k,n],
 588  f:first(args),
 589  k:second(args),
 590  n:third(args),
 591  if length(args)=3 then
 592    ZeilbergerVerboseOpt(f,k,n,verbose)
 593  else
 594    parGosperVerboseOpt(f,k,n,fourth(args),verbose)
 595  );
 596
 597
 598 /* Non verbose version of Zeilberger's algorithm */
 599 ZeilbergerVeryVerbose([args]):=
 600  block([f,k,n],
 601  f:first(args),
 602  k:second(args),
 603  n:third(args),
 604  if length(args)=3 then
 605    ZeilbergerVerboseOpt(f,k,n,very_verbose)
 606  else
 607    parGosperVerboseOpt(f,k,n,fourth(args),very_verbose)
 608  );
 609
 610
 611 /* Non verbose version of Zeilberger's algorithm */
 612 ZeilbergerExtra([args]):=
 613  block([f,k,n],
 614  f:first(args),
 615  k:second(args),
 616  n:third(args),
 617  if length(args)=3 then
 618    ZeilbergerVerboseOpt(f,k,n,extra)
 619  else
 620    parGosperVerboseOpt(f,k,n,fourth(args),extra)
 621  );
 622
 623
 624 /* Non verbose version of Zeilberger's algorithm */
 625 Zeilberger([args]):=
 626  block([f,k,n],
 627  f:first(args),
 628  k:second(args),
 629  n:third(args),
 630  if length(args)=3 then
 631    ZeilbergerVerboseOpt(f,k,n,nonverbose)
 632  else
 633    parGosperVerboseOpt(f,k,n,fourth(args),nonverbose)
 634  );
 635
 636
 637 /* Short name macro for verbose parGosper */
 638 parGosperVerbose(f,k,n,degree)::=
 639    buildq([f,k,n,degree],
 640           parGosperVerboseOpt(f,k,n,degree,verbose));
 641
 642
 643 /* Short name macro for very verbose parGosper */
 644 parGosperVeryVerbose(f,k,n,degree)::=
 645    buildq([f,k,n,degree],
 646           parGosperVerboseOpt(f,k,n,degree,very_verbose));
 647
 648
 649 /* Short name macro for the linear system version of parGosper */
 650 parGosperExtra(f,k,n,degree) ::=
 651    buildq([f,k,n,degree],
 652           parGosperVerboseOpt(f,k,n,degree,extra));
 653
 654 /* Short name macro for the summary version of parGosper */
 655 parGosperSummary(f,k,n,degree) ::=
 656    buildq([f,k,n,degree],
 657           parGosperVerboseOpt(f,k,n,degree,summary));
 658
 659
 660 /* Short name for the non-verbose version of parGosper */
 661 parGosperNonVerbose(f,k,n,degree) ::=
 662    buildq([f,k,n,degree],
 663           parGosperVerboseOpt(f,k,n,degree,nonverbose));
 664
 665 parGosper(f,k,n,order) ::=
 666    buildq([f,k,n,order],
 667           parGosperVerboseOpt(f,k,n,order,nonverbose));
 668
 669 /* Pretty print of a list */
 670 prettyPrint(list) :=
 671    block(
 672    len : length(list),
 673    for i : 1 step 1 thru len do
 674       list[i]
 675    );
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698