share/contrib/Zeilberger/makeGosperForm.mac

   1 /* Gosper form routine */
   2
   3 /* RISC Institite, Linz, Austria */
   4 /* by Fabrizio Caruso            */
   5
   6
   7 /* List select wrt a criterion */
   8 select(list,crit) :=
   9   if list = [] then
  10     []
  11   else
  12      if crit(first(list)) then
  13        cons(first(list),select(rest(list),crit))
  14      else
  15        select(rest(list),crit);
  16
  17 /* List split wrt a criterion */
  18 splitList(list,crit) :=
  19   block([trueList,falseList,i],
  20    trueList : [],
  21    falseList : [],
  22    unless list = [] do
  23       (
  24       if crit(first(list)) then
  25          trueList : endcons(first(list),trueList)
  26       else
  27          falseList : endcons(first(list),falseList),
  28       list : rest(list)
  29       ),
  30    return([trueList,falseList])
  31    );
  32
  33 /* true if a given expr has a given degree */
  34 degreeSelect(expr,var,deg) := is(degree(expr,var)=deg);
  35
  36
  37 /* integer distance */
  38 intDistance(expr1,expr2) := integerp(expr1-expr2);
  39
  40
  41 /* It splits a list into a part with a given coefficient and the rest */
  42 candSplit(list,coef,constCoef) :=
  43   block([candRes,restRes],
  44   candRes : [],
  45   restRes : [],
  46   unless (list = []) do
  47     (
  48     if (first(first(list))=coef) and integerp(second(first(list))-constCoef) then
  49        candRes : endcons(first(list),candRes)
  50     else
  51        restRes : endcons(first(list),restRes),
  52     list : rest(list)
  53     ),
  54   return([candRes, restRes])
  55   );
  56
  57
  58
  59 /* It checks whether a factored expression is linear */
  60 /* with integer coefficients */
  61 integerLinear(expr,var) :=
  62   if not(zb_operatorp(expr,"^")) then
  63     is(degree(expr,var) < 2)
  64 /*
  65 and
  66     (integerp(coeff(expr,var,1)))
  67 */
  68   else
  69     integerLinear(first(expr),var);
  70
  71 /* Separate the constant and the normalized part of an integer linear*/
  72 intLinConstSep(intLin) :=
  73   block([i,intLinRes,constRes],
  74    intLinRes : [],
  75    constRes : 1,
  76    for i : 1 thru length(intLin) do
  77       (
  78       intLinRes : endcons(first(part(intLin,i)),intLinRes),
  79       constRes : constRes * second(part(intLin,i))
  80       ),
  81    return([intLinRes, constRes])
  82   );
  83
  84
  85 constSep(intLin) :=
  86   block([i,intLinRes,constRes],
  87    intLinRes : [],
  88    constRes : 1,
  89    for i : 1 thru length(intLin) do
  90       (
  91       intLinRes : endcons(first(part(intLin,i)),intLinRes),
  92       constRes : constRes * second(part(intLin,i))
  93       ),
  94    return([intLinRes, constRes])
  95   );
  96
  97
  98 makeGosperForm(expr,k) :=
  99   makeGosperFormVerboseOpt(expr,k,nonvernose);
 100
 101 makeGosperFormVerboseOpt(expr,k,mode) :=
 102   first(makeGosperFormPrHypCheckVerboseOpt(expr,k,mode));
 103
 104 makeGosperPrHypCheck(expr,k) :=
 105   makeGosperFormPrHypCheckVerboseOpt(expr,k,nonverbose);
 106
 107 /* Check Gosper's condition */
 108 /* We are assuming that the input is factored */
 109 /* into irreducible factors */
 110 /* */
 111 /* The numerator is q(k) and the denominator is r(k) */
 112 makeGosperFormPrHypCheckVerboseOpt(expr,k,mode) :=
 113   block(
 114   /* Local variables */
 115   [num,den,numList,denList,lengthNumList,lengthDenList,
 116    test,
 117    i__,j,h,minExp,pulledPoly,pPoly,curExp,compPoly,formSign,
 118    p,q,r,pList,qList,rList,
 119    negNumFlag,negDenFlag,
 120    intLinNum,nonIntLinNum,intLinDen,nonIntLinDen,
 121    sepNumList,sepDenList,
 122    constIntLinNum, constIntLinDen, constSep,
 123    normIntLinNum, normNonIntLinNum, normIntLinDen, normNonIntLinDen,
 124    numItem, denItem,
 125    numExp,denExp,diff,candListSplit,item,hyp_flag, candNumList, candDenList, intLinNumConst, intLinDenConst, nonIntLinNumConst, nonIntLinDenConst],
 126
 127   hyp_flag : true,
 128   if(expr=1) then
 129     return([[1,1,1],true]),
 130
 131 /*
 132     num:factor(num(expr)),
 133     den:factor(denom(expr)),
 134 */
 135     num : num(expr),
 136     den : denom(expr),
 137
 138
 139   if mode>= verbose then
 140     (
 141     print(num),
 142     print(den)
 143     ),
 144
 145 negNumFlag : zb_operatorp(num,"-"),
 146 negDenFlag : zb_operatorp(den,"-"),
 147
 148 if (negNumFlag) then
 149   num : -1 * num,
 150 if (negDenFlag) then
 151   den : -1 * den,
 152
 153 if (negNumFlag and not(negDenFlag))or(not(negNumFlag) and negDenFlag) then
 154    (
 155    formSign : -1
 156    )
 157 else
 158    formSign : 1,
 159
 160 pPoly : 1,
 161
 162
 163 sepNumList:intLinSep(num,k),
 164 sepDenList:intLinSep(den,k),
 165
 166 if mode >= verbose then
 167 (
 168 print("sepNumList : ", sepNumList),
 169 print("sepDenList : ", sepDenList)
 170 ),
 171
 172 nonIntLinNum : first(sepNumList),
 173 intLinNum : second(sepNumList),
 174 nonIntLinDen : first(sepDenList),
 175 intLinDen : second(sepDenList),
 176
 177 p:1,
 178 q:1,
 179 r:1,
 180 pList:[],
 181 qList:[],
 182 rList:[],
 183
 184 if numList # [] or denList # [] then
 185   hyp_flag : false,
 186
 187 /* Non-integer linear computation */
 188
 189 numList : normList(nonIntLinNum,k),
 190 denList : normList(nonIntLinDen,k),
 191
 192 if mode>= verbose then
 193 (
 194 print("numList : ", numList),
 195 print("denList : ", denList)
 196 ),
 197
 198 constSep : constSep(numList),
 199 normNonIntLinNum : constSep[1],
 200 nonIntLinNumConst : constSep[2],
 201 if mode>= verbose then
 202 (
 203 print("normIntLinNum : ", normNonIntLinNum),
 204 print("intLinNumConst : ", nonIntLinNumConst)
 205 ),
 206
 207 constSep : constSep(denList),
 208 normNonIntLinDen : constSep[1],
 209 nonIntLinDenConst : constSep[2],
 210 if mode>= verbose then
 211 (
 212 print("normIntLinDen : ", normNonIntLinDen),
 213 print("nonIntLinDenConst : ", nonIntLinDenConst)
 214 ),
 215
 216 q : q * nonIntLinNumConst,
 217 r : r * nonIntLinDenConst,
 218
 219 numList : normNonIntLinNum,
 220 denList : normNonIntLinDen,
 221
 222 /* entering outer loop */
 223
 224 if mode >= verbose then
 225   print("denList : ", denList),
 226 unless ((numList = []) or (denList = [])) do
 227   (
 228   denItem : first(denList),
 229   if mode >= verbose then
 230      print("denItem : ", denItem),
 231   candListSplit : splitList(numList,
 232                             lambda([x],degreeSelect(x,k,degree(first(denList),k)))
 233                            ),
 234 /*
 235 candListSplit : splitList(numList,
 236                             lambda([x],degreeSelect(x,k,degree(first(denItem),k)))
 237                            ),
 238 */
 239
 240   candNumList : first(candListSplit),
 241
 242 numList : second(candListSplit),
 243
 244 candListSplit : splitList(denList,
 245                             lambda([x],degreeSelect(x,k,degree(first(denList),k)))
 246                            ),
 247 /*
 248 candListSplit : splitList(denList,
 249                             lambda([x],degreeSelect(x,k,degree(first(denItem),k)))
 250                            ),
 251 */
 252
 253
 254 candDenList : first(candListSplit),
 255 denList : second(candListSplit),
 256
 257
 258
 259 /*  inner loop */
 260   unless candDenList = [] do
 261      (
 262      denItem : first(candDenList),
 263      tempList : candNumList,
 264      candNumList : [],
 265      unless tempList = [] do
 266         (
 267
 268
 269      shiftRes : resultant(norm2polyPower(first(tempList)),
 270                           expand(subst(k+j,k,norm2polyPower(denItem))),k),
 271      if mode >= very_verbose then
 272         print("shiftRes : ", shiftRes),
 273
 274      factList : poly2list(factor(shiftRes)),
 275
 276      intFactList : first(splitList(factList, lambda([y],integerLinear(y,j)and degree(y,j)>0))),
 277
 278      solList : makelist(second(first(solve(part(intFactList,i),j))),i,1,length(intFactList)),
 279
 280      intSolList : first(splitList(solList,lambda([y],integerp(y) and (y>=0)))),
 281
 282      if intSolList # [] then
 283          (
 284          minSol : apply(min,intSolList),
 285
 286          numExp : second(first(tempList)),
 287          denExp : second(denItem),
 288          minExp : min(numExp,denExp),
 289
 290          if numExp = denExp then
 291              (
 292              p : p * product(norm2polyPower(subst(k+j,k,denItem)),j,0,minSol-1),
 293
 294              denItem : [1,1],
 295              tempList : []
 296
 297              )
 298          else
 299              (
 300              item : [first(denItem),minExp],
 301
 302              p : p * product(norm2polyPower(subst(k+j,k,item)),j,0,minSol-1),
 303
 304              if numExp > denExp then
 305                 (
 306                 candNumList : endcons([first(first(tempList)),numExp-minExp],candNumList),
 307                 denItem : [1,1],
 308                 tempList : []
 309                 )
 310              else
 311                 (
 312                 denItem : [first(denItem), denExp-minExp],
 313                 tempList : rest(tempList)
 314                 )
 315              )
 316          )
 317      else
 318          (
 319          if mode>= verbose then
 320             print("no positive integer solutions"),
 321          candNumList : endcons(first(tempList),candNumList),
 322          tempList : rest(tempList)
 323          )
 324
 325         ),
 326      r : r * norm2polyPower(denItem),
 327      candDenList : rest(candDenList)
 328
 329      ),
 330
 331 /* end of inner loop */
 332 q : q * product(norm2polyPower(part(candNumList,i__)),i__,1,length(candNumList))
 333
 334   ),
 335 /* end of outer loop */
 336 if mode>= verbose then
 337 (
 338 print("numList : ", numList),
 339 print("denList : ", denList)
 340 ),
 341 if numList # [] then
 342   q : q * product(norm2polyPower(part(numList,i__)),i__,1,length(numList)),
 343 if denList # [] then
 344   r : r * product(norm2polyPower(part(denList,i__)),i__,1,length(denList)),
 345
 346 if mode>= verbose then
 347 (
 348 print("q : ", q),
 349 print("r : ", r)
 350 ),
 351
 352
 353 /* Integer linear computation */
 354
 355 numList : intLinNormList(intLinNum,k),
 356
 357 numList : sort(numList),
 358
 359 denList : intLinNormList(intLinDen,k),
 360
 361 denList : sort(denList),
 362
 363 constSep : intLinConstSep(numList),
 364 normIntLinNum : constSep[1],
 365 intLinNumConst : constSep[2],
 366 if mode>= verbose then
 367 (
 368 print("normIntLinNum : ", normIntLinNum),
 369 print("intLinNumConst : ", intLinNumConst)
 370 ),
 371
 372 constSep : intLinConstSep(denList),
 373 normIntLinDen : constSep[1],
 374 intLinDenConst : constSep[2],
 375 if mode>= verbose then
 376 (
 377 print("normIntLinDen : ", normIntLinDen),
 378 print("intLinDenConst : ", intLinDenConst)
 379 ),
 380
 381 numList : normIntLinNum,
 382 denList : normIntLinDen,
 383
 384   lengthNumList:length(numList),
 385   lengthDenList:length(denList),
 386
 387 unless ((numList = []) or (denList = [])) do
 388   (
 389   denItem : first(denList),
 390
 391 /* This is an adaptation of Schorn-Riese's Mathematica code */
 392
 393 if mode>= verbose then
 394 (
 395 print("denItem : ", denItem),
 396 print("numList : ", numList),
 397 print("denList : ", denList)
 398 ),
 399
 400
 401 candListSplit : splitList(numList,
 402                 lambda([x],intDistance(first(x),first(denItem)))),
 403
 404 if mode>= verbose then
 405   print("num) candListSplit : ", candListSplit),
 406
 407
 408 candNumList : first(candListSplit),
 409
 410 numList : second(candListSplit),
 411 candListSplit : splitList(denList,
 412                  lambda([x],intDistance(first(x),first(denItem)))),
 413
 414 if mode>= verbose then
 415   print("den) candListSplit : ", candListSplit),
 416
 417
 418 candDenList : first(candListSplit),
 419
 420 denList : second(candListSplit),
 421
 422
 423 unless ((candNumList = []) or (candDenList = [])) do
 424     (
 425
 426     diff : first(first(candNumList))-first(first(candDenList)),
 427     if mode>= verbose then
 428       (
 429       print("diff : ", diff),
 430       print("inner cycle) candNumList : ", candNumList),
 431       print("inner cycle) candDenList : ", candDenList)
 432       ),
 433     if not(integerp(diff)) or (integerp(diff) and (diff<0)) then
 434        (
 435        /* q must be updated here */
 436        /* Gosper-condition satisfied */
 437
 438        q : q * norm2polyPower(first(candNumList)),
 439        candNumList : rest(candNumList)
 440        )
 441     else
 442        (
 443        /* Gosper-condition NOT satisfied */
 444        /* p must be updated here */
 445        numExp:second(first(candNumList)),
 446        denExp:second(first(candDenList)),
 447        minExp:min(numExp,denExp),
 448
 449        p : p * product(subst(k+i__,k,norm2intLin(first(candDenList),k)^minExp),i__,0,diff-1),
 450        if mode>= verbose then
 451           print("p computed : ", p),
 452        if(numExp>denExp) then
 453           (
 454           item : [first(first(candNumList)),numExp-minExp],
 455
 456           candNumList : cons(item,rest(candNumList)),
 457
 458           candDenList : rest(candDenList)
 459           )
 460        else
 461           if(numExp<denExp) then
 462              (
 463              item : [first(first(candDenList)),denExp-minExp],
 464              candNumList : rest(candNumList),
 465              candDenList : cons(item,rest(candDenList))
 466              )
 467           else
 468             (
 469             candNumList : rest(candNumList),
 470             candDenList : rest(candDenList)
 471             )
 472
 473        )
 474     ), /* end inner unless cycle */
 475
 476 /* out of inner cycle */
 477
 478 if (candNumList # []) then
 479  q : q * product(norm2polyPower(part(candNumList,i__)),i__,1,length(candNumList)),
 480
 481
 482 if(candDenList # []) then
 483  r : r * product(norm2polyPower(part(candDenList,i__)),i__,1,length(candDenList))
 484
 485  ), /* end unless */
 486
 487 if (numList # []) then
 488   q : q*product(norm2polyPower(part(numList,i__)),i__,1,length(numList)),
 489
 490 if (denList # []) then
 491   r : r*product(norm2polyPower(part(denList,i__)),i__,1,length(denList)),
 492
 493
 494 r : subst(k-1,k,r),
 495
 496 return([[p,q*formSign*intLinNumConst,r*intLinDenConst],hyp_flag])
 497
 498
 499 ); /* end block */
 500
 501
 502
 503
 504 makeGosperForm(expr,var) ::=
 505   buildq([expr,var],makeGosperFormVerboseOpt(expr,var,nonverbose));
 506
 507 makeGosperFormVerbose(expr,var) ::=
 508   buildq([expr,var],makeGosperFormVerboseOpt(expr,var,verbose));
 509
 510
 511
 512
 513 /* Macro that extracts p out of the Gosper form */
 514 takeP(GosperForm) ::=
 515   buildq([GosperForm],first(GosperForm));
 516
 517
 518 /* Macro that extracts q out of the Gosper form */
 519 takeQ(GosperForm) ::=
 520   buildq([GosperForm],second(GosperForm));
 521
 522
 523 /* Macro that extracts r out of the Gosper form */
 524 takeR(GosperForm) ::=
 525   buildq([GosperForm],third(GosperForm));
 526
 527
 528
 529
 530
 531