share/contrib/Zeilberger/shiftQuotient.mac

   1 /* Routines for computing:  */
   2 /* a) application of a lin. op. on a hyp. seq. */
   3 /* b) shift quotient of a hyp. seq. */
   4
   5 /* RISC Institute, Linz, Austria */
   6 /* by Fabrizio Caruso            */
   7
   8
   9 dependent(expr,var) :=
  10   block([res,count],
  11     if(expr=var) then
  12       return(true)
  13     else
  14       if(atom(expr)) then
  15         return(false)
  16       else
  17         (
  18         res:false,
  19
  20         for count:1 unless (res=true)or(count>length(expr)) do
  21            (
  22            res : dependent(part(expr,count),var)
  23
  24            ),
  25         return(res)
  26         )
  27 );
  28
  29
  30 poly2list(pol) :=
  31    block(
  32    if atom(pol) or not(zb_operatorp(pol,"*")) then
  33      return([pol])
  34    else
  35      return(cons(first(pol),poly2list(pol/first(pol))))
  36    );
  37
  38 extractConstant(polyList,var) :=
  39    block([resConst, resDep,i],
  40    resDep : 1,
  41    resConst : 1,
  42    for i : 1 thru length(polyList) do
  43       if (atom(polyList[i]) and polyList[i]=var) or
  44          (not(atom(polyList[i])) and dependent(polyList[i],var)) then
  45          resDep : resDep * polyList[i]
  46       else
  47          resConst : resConst * polyList[i],
  48    return([resDep,resConst])
  49    );
  50
  51
  52 /* It computes the rat. factor out of the appl. of a lin. op. to a hyp. seq. */
  53 niceForm(hyp,var,parName,ord,sumVar) :=
  54     block(
  55     [shQuo,numConst,denConst,num,den,res,i],
  56
  57     res:parName[ord],
  58     shQuo : shiftQuoHypCheck(hyp,var),
  59     if second(shQuo) = NO_HYP then
  60       return([])
  61     else shQuo : first(shQuo),
  62
  63     numConst : extractConstant(poly2list(factor(num(shQuo))),sumVar),
  64     denConst : extractConstant(poly2list(factor(denom(shQuo))),sumVar),
  65
  66     num : numConst[1],
  67     den : denConst[1],
  68
  69     shQuo : num/den,
  70     for i : ord step -1 thru 1 do
  71       res : xthru(parName[i-1] + factor(shiftFactPoly(shQuo,var,i-1))*res),
  72
  73     return([res,numConst[2],denConst[2]])
  74 );
  75
  76
  77 /* It computes the rat. factor out of the appl. of a lin. op. to a hyp. seq. */
  78 niceFormDB(hyp,var,parName,ord) :=
  79     block(
  80     [shQuo,num,den,res,i],
  81
  82     res:parName[ord],
  83     shQuo : shiftQuo(hyp,var),
  84     print("Shift quotient computed!"),
  85     print("Order : ", ord),
  86     for i : ord step -1 thru 1 do
  87       (
  88       res : xthru(parName[i-1] + shiftFactPoly(shQuo,n,i-1)*res),
  89       print("i : ", i)
  90       ),
  91     return(res)
  92     );
  93
  94
  95 /*
  96 removeBinomial(expr) :=
  97   if atom(expr) then
  98     expr
  99   else
 100     if zb_op(expr) = binomial then
 101       first(expr)!/second(expr)!/(first(expr)-second(expr))!
 102     else
 103       apply(zb_op(expr),makelist(removeBinomial(part(expr,i)),i,1,length(expr)));
 104 */
 105
 106 removeBinomial(expr) :=
 107   block([],
 108
 109   if mapatom(expr) then
 110     return(expr)
 111   else
 112     if zb_op(expr) = binomial then
 113       return(first(expr)!/second(expr)!/(first(expr)-second(expr))!)
 114     else
 115       if zb_op(expr) = "-" then
 116         return(-removeBinomial(first(expr)))
 117       else
 118
 119         map(removeBinomial, expr)
 120   );
 121
 122 shiftFactPoly(expr,k,j) :=
 123   if atom(expr) then
 124     subst(k+j,k,expr)
 125   else
 126   if zb_op(expr) = "/" then
 127      shiftFactPoly(first(expr),k,j)*shiftFactPoly(second(expr),k,j)^(-1)
 128   else
 129   if zb_op(expr) = "*" then
 130      product(shiftFactPoly(part(expr,i__),k,j),i__,1,length(expr))
 131   else
 132   if zb_op(expr) = "^" then
 133      if integerp(second(expr)) then
 134          shiftFactPoly(first(expr),k,j)^second(expr)
 135      else
 136         0
 137   else
 138     expand(subst(k+j,k,expr));
 139
 140
 141 shiftQuo(expr,k) :=
 142   block([sq_res],
 143     sq_res : shiftQuoHypCheck(expr,k),
 144     if second(sq_res) = HYP then
 145       return(first(sq_res))
 146     else
 147       return(subst(k+1,k,first(sq_res))/first(sq_res))
 148     );
 149
 150
 151 shiftQuoOnlyHyp(expr,k) :=
 152   block([sq_res],
 153   sq_res : shiftQuoHypCheck(expr,k),
 154   if second(sq_res) = HYP then
 155     return(first(sq_res))
 156   else
 157     return([])
 158   );
 159
 160
 161 isPolynomial(expr,k) :=
 162   if freeof(k,expr) or expr=k or expr=-k or constantp(expr) then
 163      true
 164   else
 165     if zb_op(expr) = "^" then
 166       ( freeof(k,second(expr)) and isPolynomial(first(expr),k) )
 167     else
 168       if zb_op(expr) = "-" then
 169         isPolynomial(first(expr),k)
 170       else
 171         if zb_op(expr) = "*" or zb_op(expr)= "+" or zb_op(expr) = "-" then
 172            apply("and", makelist(isPolynomial(part(expr,i),k),i,1,length(expr)))
 173         else
 174            false;
 175
 176
 177
 178 rationalp(expr,k) :=
 179   if not(isPolynomial(expr,k)) and
 180      not(zb_op(expr) = "/" and
 181                  isPolynomial(expand(num(expr)),k) and
 182                  isPolynomial(expand(denom(expr)), k)) then
 183      (
 184      if zb_op(expr) = "-" then
 185        rationalp(first(expr),k)
 186      else
 187        false
 188      )
 189   else
 190      true;
 191
 192 shiftQuoHypCheck(expr,k) :=
 193   block([xthru_expr,sq_res],
 194     xthru_expr : xthru(removeBinomial(expr)),
 195     sq_res : shiftQuoHypCheckAux(xthru_expr,k,HYP),
 196
 197     if rationalp(sq_res[1],k) then
 198       return(sq_res)
 199     else
 200       return([sq_res[1],NO_HYP])
 201
 202 );
 203
 204
 205 shiftQuoHypCheckAux(expr,k,hyp_flag) :=
 206 block([sq_num,sq_den,sq_base,sq_exp],
 207   if hyp_flag = NO_HYP then
 208     return([expr,NO_HYP])
 209   else
 210   if freeof(k,expr) then
 211      return([1,HYP])
 212   else
 213      if expr = k then
 214        return([(k+1)/k,HYP])
 215      else
 216      if zb_op(expr) = "*" then
 217        return(product(shiftQuoHypCheckAux(part(expr,i__),k,hyp_flag),
 218               i__,1,length(expr)))
 219      else
 220      if zb_op(expr) = "/" then
 221        (
 222        sq_num : shiftQuoHypCheckAux(first(expr),k,hyp_flag),
 223        sq_den : shiftQuoHypCheckAux(second(expr),k,hyp_flag),
 224        return([first(sq_num)/first(sq_den),second(sq_num)*second(sq_den)])
 225        )
 226      else
 227      if zb_op(expr) = "^" then
 228        (
 229        if (freeof(k, second(expr))) then
 230          (
 231          sq_base : shiftQuoHypCheckAux(first(expr), k, hyp_flag),
 232          if sq_base[2]=NO_HYP then return([expr, NO_HYP])
 233          else return([sq_base[1]^second(expr), hyp_flag])
 234          )
 235        else if (freeof(k, first(expr))) then
 236          (
 237          sq_exp: bothcoef(expand(second(expr)), k),
 238          if not(freeof(k, second(sq_exp))) then return([expr, NO_HYP])
 239          else return([first(expr)^first(sq_exp), hyp_flag])
 240          )
 241        else return([expr, NO_HYP])
 242        )
 243      else
 244      if zb_op(expr) = "!" then
 245         if not(integerp(leadCoeff(first(expr),k))) then
 246            return([expr,NO_HYP])
 247         else
 248           if leadCoeff(first(expr),k)>0 then
 249             return(
 250                 [product(factor(first(expr)+i__),i__,1,leadCoeff(first(expr),k)),
 251                 hyp_flag])
 252           else
 253              [
 254              1/product(factor(first(expr)-i__+1),
 255                      i__,1,-leadCoeff(first(expr),k)),hyp_flag]
 256 /*
 257      else
 258      if zb_op(expr) = binomial then
 259         (
 260         sq_num : shiftQuoHypCheckAux(factorial(first(expr)),k,hyp_flag),
 261         sq_den : shiftQuoHypCheckAux(factorial(first(expr)-second(expr)),
 262                                   k,hyp_flag)*
 263                  shiftQuoHypCheckAux(factorial(second(expr)),
 264                                   k,hyp_flag),
 265         if second(sq_num) = HYP and second(sq_den)=HYP then
 266           return([first(sq_num)/first(sq_den),HYP])
 267         else
 268           return([expr,NO_HYP])
 269
 270         )
 271 */
 272      else
 273        if zb_op(expr) = "-" then block(
 274          return(shiftQuoHypCheckAux(-expr, k, hyp_flag))
 275        )
 276        else
 277        if zb_op(expr) = "+" then
 278          return([shiftFactPoly(expr,k,1)/expr,hyp_flag])
 279        else
 280          (
 281          if warnings then print("Unknown operator : ", zb_op(expr)),
 282          return([expr,NO_HYP])
 283          )
 284 );
 285
 286
 287
 288
 289
 290
 291
 292