share/contrib/Zeilberger/GosperEq.mac

   1
   2 /* Gosper's equation-related routines */
   3
   4 /* RISC Institite, Linz, Austria */
   5 /* by Fabrizio Caruso            */
   6
   7
   8
   9
  10
  11 /* It creates the Gosper's polynomial out of the following parameters:*/
  12 /* polynomials p,q,r + */
  13 /* Ansatz degree deg of the solution + */
  14 /* indeterminate + */
  15 /* unknown coefficients name */
  16 AnsatzDegreeVerboseOpt(p,q,r,indet,mode) :=
  17   block(
  18   [qp,qm,degqp,degqm,degp,qpCoeff,qmCoeff,alpha],
  19
  20   /* Initialization */
  21   qp : expand(q + r),
  22   qm : expand(q - r),
  23   degqp : degree(expand(qp),indet),
  24   degqm : degree(expand(qm),indet),
  25   degp : degree(expand(p),indet),
  26
  27   if mode >= verbose then
  28     (
  29     print("(Ansatz) Degree of qplus:" , degqp),
  30     print("(Ansatz) Degree of qminus: " , degqm),
  31     print("(Ansatz) Degree of p: " , degp)
  32     ),
  33
  34   /* Case distinction */
  35   if degqm >= degqp then
  36      (
  37      if mode >= verbose then
  38        (
  39        print("(Ansatz degree) Ordinary case (degqm>=degqp)"),
  40        print("(Ansatz degree) Degree : ", degp-degqm)
  41        ),
  42      return(degp - degqm)
  43      )
  44   else
  45      (
  46      if mode >= verbose then
  47        print("(Ansatz degree) Non-ordinary case (degqm<degqp)"),
  48      qpCoeff : leadCoeff(qp,indet),
  49      if mode >= verbose then
  50        print("(Ansatz degree) Leading coeff. of qplus: ", qpCoeff),
  51      qmCoeff : leadCoeff(qm,indet),
  52      if mode >= verbose then
  53        print("(Ansatz degree) Leading coeff. of qminus: ", qmCoeff),
  54      alpha : -2*(qmCoeff/qpCoeff),
  55      if mode >= verbose then
  56        print("(Ansatz degree) alpha :", alpha),
  57      if integerp(alpha) then
  58         (
  59         if mode >= verbose then
  60           print("(Ansatz degree) Nasty case (int(-2*qmC/qpC))"),
  61         return(max(alpha,degp-degqp+1,0)) /* 0??? */
  62         )
  63      else
  64         (
  65         if mode >= verbose then
  66           print("(Ansatz degree) Easy case (not(int(-2*qmC/qpC)))"),
  67         return(max(degp-degqp+1,0)) /* 0??? */
  68         )
  69      )
  70    );
  71
  72
  73
  74
  75 /* p,q,r: Gosper's form polynomials */
  76 /* indet: indeterminate */
  77 /* coef: prefix of the coefficients */
  78 /* deg: desidered degree of the Gosper's polynomial */
  79 GosperPoly(p,q,r,indet,coef,deg) :=
  80    expand(q*subst(k+1,k,poly(indet,coef,deg))-r*poly(indet,coef,deg)-p);
  81
  82
  83
  84
  85 /* It solves the Gosper's equation (Gosper's polynomial=0) */
  86 /* and yields the coefficients of the solution in a list */
  87 /* */
  88 /* We are assuming that Gpoly is fully expanded */
  89
  90
  91 /* It solves the Gosper's equation (Gosper's polynomial=0) */
  92 /* and yields the coefficients of the solution in a list */
  93 /* */
  94 /* We are assuming that Gpoly is fully expanded */
  95
  96 /* Gpoly: Gosper polynomial */
  97 /* unkPar: */
  98 /* numUnkPar: */
  99 /* indet : */
 100 SolveSysVerboseOpt(Gpoly,unkPar,numUnkPar,indet,mode) :=
 101    block([dim,base,baseElem,sol],
 102
 103    if mode >= verbose then
 104      (
 105      print("(SolveSys) Poly2List: ", poly2List(Gpoly,indet)),
 106      print("(SolveSys) Variables: ",
 107             makelist(unkPar[i],i,0,numUnkPar-1))
 108      ),
 109
 110   sol : linear_solver(poly2List(Gpoly,indet),
 111                       makelist(unkPar[i],i,0,numUnkPar-1)),
 112
 113   dim : length(ind_vars(linear_solver)),
 114   if mode >= verbose then
 115     (
 116     print("mode : ", mode),
 117     print("dim : ", dim),
 118     print("debug : ", ind_vars(linear_solver)),
 119     print("sol : ", sol)
 120     ),
 121
 122   if dim>0 and not(sol=[]) then
 123      (
 124      base : [],
 125      for i:1 thru dim do
 126        (
 127        baseElem : sol,
 128        for j:1 thru dim do
 129           if j = i then
 130              baseElem : subst(0,ind_vars(linear_solver)[j],baseElem)
 131           else
 132              baseElem : subst(1,ind_vars(linear_solver)[j],baseElem),
 133        base : endcons(baseElem,base)
 134        ),
 135      return(base)
 136      )
 137   else
 138      return([sol])
 139    );
 140
 141
 142 /* It solves the Gosper's equation (Gosper's polynomial=0) */
 143 /* and yields the coefficients of the solution in a list */
 144 /* */
 145 /* We are assuming that Gpoly is fully expanded */
 146
 147 /* Gpoly: Gosper polynomial */
 148 /* unkPar: */
 149 /* numUnkPar: */
 150 /* indet : */
 151 SolveZSysVerboseOpt(Gpoly,unkPar1,numUnkPar1,unkPar2,numUnkPar2,indet,
 152                     solver,mode) :=
 153    block([sol,dim,base,i,j,baseElem,sys],
 154
 155    if mode>=verbose then
 156      (
 157      print("inside : ", mode),
 158      print("(SolveZSys) Unknowns: ", unkPar1, unkPar2),
 159      print("(SolveZSys) Number of unknowns: ", numUnkPar1, numUnkPar2),
 160      print("(SolveZSys) Gpoly : ", Gpoly),
 161      print("(SolveZSys) Poly2List: ", poly2List(Gpoly,indet)),
 162      print("(SolveZSys) Variables: ",
 163          append(
 164          makelist(unkPar1[i],i,0,numUnkPar1-1),
 165          makelist(unkPar2[i],i,0,numUnkPar2-1))),
 166      print("Gpoly : ", Gpoly)
 167      ),
 168    sys : poly2List(Gpoly,indet),
 169
 170    if mode=verbose then
 171      print("sys : ", sys),
 172
 173    sys : makelist(clearGCD(part(sys,i),unkPar1,numUnkPar1,
 174                                        unkPar2,numUnkPar2),i,1,length(sys)),
 175
 176    sol : solver(poly2List(Gpoly,indet),
 177          append(
 178          makelist(unkPar1[i],i,0,numUnkPar1-1),
 179          makelist(unkPar2[i],i,0,numUnkPar2-1))),
 180    if mode=verbose then
 181      (
 182      print("sys : ", sys),
 183      print("sol : ", sol)
 184      ),
 185    dim : length(ind_vars(solver)),
 186    base : [],
 187    for i:1 thru dim do
 188       (
 189       baseElem : sol,
 190       for j:1 thru dim do
 191          if j = i then
 192             baseElem : subst(1,ind_vars(solver)[j],baseElem)
 193          else
 194             baseElem : subst(0,ind_vars(solver)[j],baseElem),
 195       base : endcons(baseElem,base)
 196       ),
 197    return(base)
 198         );
 199
 200 checkSol(sys,sol) :=
 201   block([i,res],
 202 res : sys,
 203 for i:1 thru length(sol) do
 204    res : subst(second(sol[i]),first(sol[i]),res),
 205 return(factor(res))
 206 );
 207
 208
 209 /* It solves the Gosper's equation (Gosper's polynomial=0) */
 210 /* and yields the coefficients of the solution in a list */
 211 /* */
 212 /* We are assuming that Gpoly is fully expanded */
 213
 214 /* Gpoly: Gosper polynomial */
 215 /* unkPar: */
 216 /* numUnkPar: */
 217 /* indet : */
 218 SolveZSysVerboseSys(Gpoly,unkPar1,numUnkPar1,unkPar2,numUnkPar2,indet) :=
 219    block(
 220    print("(SolveZSys) Unknowns: ", unkPar1, unkPar2),
 221    print("(SolveZSys) Number of unknowns: ", numUnkPar1, numUnkPar2),
 222    print("(SolveZSys) Poly2List: ", poly2List(Gpoly,indet)),
 223    print("First equation : ", part(poly2List(Gpoly,indet),1)),
 224    print("Coeff in _z[0] of degree 0 : ", coeff(part(poly2List(Gpoly,indet),1),n,0)),
 225    print("(SolveZSys) Variables: ",
 226          append(
 227          makelist(unkPar1[i],i,0,numUnkPar2-1),
 228          makelist(unkPar2[i],i,0,numUnkPar1-1))),
 229    print(""),
 230    print("(SolveZSys) LINEAR SYSTEM OF EQUATIONS"),
 231    showSys(Gpoly,indet),
 232
 233    [linear_solver(poly2List(Gpoly,indet),
 234          append(
 235          makelist(unkPar1[i],i,0,numUnkPar2-1),
 236          makelist(unkPar2[i],i,0,numUnkPar1-1)))]
 237         );
 238
 239
 240
 241 /* It shows the corresponding linear system of a Gosper-polynomial */
 242 showSys(Gpoly,indet) :=
 243    printSys(poly2List(Gpoly,indet),1);
 244
 245
 246 /* Pretty print of am homogeneus system of linear equations */
 247 printSys(list,num) :=
 248    block(
 249
 250    print(""),
 251    if length(list) = 1 then
 252       print(first(list)," = 0.    (",num,")" )
 253    else
 254       (
 255       print(""),
 256       print(first(list)," = 0,    (",num,")" ),
 257       printSys(rest(list),num+1)
 258       )
 259    );
 260
 261
 262 /* Pretty print of the degrees of the solutions of a system of linear equations */
 263 /* Separates the Gosper and Zeilberger Ansatz coefficients */
 264 getDegrees(Gres,GLength,ZLength,n) :=
 265   block(
 266   [GDegList,ZDegList],
 267
 268   minListLength : min(GLength,ZLength),
 269   GDegList:[],
 270   ZDegList:[],
 271
 272   for i : 1 thru minListLength*2 do
 273      (
 274      if oddp(i) then
 275         (
 276         print(lhs(Gres[i]), "=", factor(xthru(rhs(Gres[i]))), "  deg: ",
 277               degree(expand(factor(xthru(rhs(Gres[i])))),n)),
 278         print(""),
 279         GDegList : endcons(expand(rhs(Gres[i])),GDegList)
 280         ),
 281
 282      if evenp(i) then
 283         (
 284         print(lhs(Gres[i]), "=", factor(xthru(rhs(Gres[i]))), "  deg: ",
 285               degree(expand(factor(xthru(rhs(Gres[i])))),n)),
 286         print(""),
 287         ZDegList : endcons(expand(rhs(Gres[i])),ZDegList)
 288         )
 289      ),
 290   if GLength > ZLength then
 291      for i: minListLength*2+1 thru GLength+ZLength do
 292          (
 293          print(lhs(Gres[i]), "=", factor(xthru(rhs(Gres[i]))), "  deg: ",
 294                degree(expand(factor(xthru((rhs(Gres[i]))))),n)),
 295          print(""),
 296          GDegList : endcons(expand(rhs(Gres[i])),GDegList)
 297          ),
 298   if ZLength > GLength then
 299      for i: minListLength*2+1 thru GLength+ZLength do
 300          (
 301          print(lhs(Gres[i]), "=", factor(xthru(rhs(Gres[i]))), "  deg: ",
 302                degree(expand(factor(xthru(rhs(Gres[i])))),n)),
 303          print(""),
 304          ZDegList : endcons(expand(rhs(Gres[i])),ZDegList)
 305          ),
 306
 307   [GDegList,ZDegList]
 308   );
 309
 310
 311
 312 /* It solves the Gosper's equation (Gosper's polynomial=0) */
 313 /* and yields the coefficients of the solution in a list */
 314 /* */
 315 /* We are assuming that Gpoly is fully expanded */
 316
 317 /* Gpoly: Gosper polynomial */
 318 /* unkPar: */
 319 /* numUnkPar: */
 320 /* indet : */
 321
 322
 323 SolveZSys(Gpoly,unkPar1,numUnkPar1,unkPar2,numUnkPar2,indet) :=
 324    block([sol,dim,base,i,j,baseElem,pList,sys,clearedSys],
 325    sys : poly2List(Gpoly,indet),
 326
 327
 328    sys : makelist(clearGCD(part(sys,i),unkPar1,numUnkPar1,unkPar2,numUnkPar2),i,1,length(sys)),
 329
 330    sol : linear_solver(sys,
 331          append(
 332          makelist(unkPar1[i],i,0,numUnkPar1-1),
 333          makelist(unkPar2[i],i,0,numUnkPar2-1))),
 334
 335    dim : length(ind_vars(linear_solver)),
 336    base : [],
 337    for i:1 thru dim do
 338       (
 339       baseElem : sol,
 340       for j:1 thru dim do
 341          if j = i then
 342             baseElem : subst(1,ind_vars(linear_solver)[j],baseElem)
 343          else
 344             baseElem : subst(0,ind_vars(linear_solver)[j],baseElem),
 345       base : endcons(baseElem,base)
 346       ),
 347    return(base)
 348     );
 349
 350 checkSol(sys,sol) :=
 351   block([i,res],
 352 res : sys,
 353 for i:1 thru length(sol) do
 354    res : subst(second(sol[i]),first(sol[i]),res),
 355 return(factor(res))
 356 );
 357
 358
 359
 360
 361
 362
 363
 364
 365