share/contrib/Zeilberger/algUtil.mac

   1 /* Some algebraic/combinatorial auxiliary routines */
   2
   3 /* RISC Institute, Linz, Austria */
   4 /* by Fabrizio Caruso            */
   5
   6
   7 /* Least common multiple */
   8 lcm(a,b) := a*b / gcd(a,b);
   9
  10
  11 /* Remove common gcd */
  12 clearGCD(pol,unkPar1,numUnkPar1,unkPar2,numUnkPar2) :=
  13   block([pList,i,cmfact,res],
  14
  15   pList : [],
  16   for i:1 thru numUnkPar1 do
  17     pList : endcons(coeff(pol,unkPar1[i-1],1),pList),
  18   for i:1 thru numUnkPar2 do
  19     pList : endcons(coeff(pol,unkPar2[i-1],1),pList),
  20   cmfact : first(pList),
  21   for i : 2 thru numUnkPar1 + numUnkPar2 - 1 do
  22     cmfact : gcd(cmfact,part(pList,i)),
  23
  24
  25   if (cmfact # 1) and (cmfact # 0) then
  26    (
  27    res : 0,
  28    for i : 1 thru numUnkPar1 do
  29      res : res + factor(part(pList,i)/cmfact) * unkPar1[i-1],
  30    for i : 1 thru numUnkPar2 do
  31      res : res + factor(part(pList,numUnkPar1+i)/cmfact) * unkPar2[i-1],
  32
  33    return(res)
  34    )
  35   else
  36    return(pol)
  37   );
  38
  39
  40 /* Make a polynomial primitive */
  41 primitive(pol,var) :=
  42   block([polyList,coeGcd,i,res],
  43     pol : expand(pol),
  44     polyList : makelist(coeff(pol,var,i),i,0,degree(pol,var)),
  45     coeGcd : part(polyList,1),
  46     for i : 2 thru length(polyList) do
  47       coeGcd : gcd(coeGcd,part(polyList,i)),
  48      res : 0,
  49     for i : 1 thru length(polyList) do
  50       res : res + factor(part(polyList,i)/coeGcd) * var^(i-1),
  51     return(res)
  52     );
  53
  54
  55 /* Make a list out of a product of powers of integers disregarding powers */
  56 factorsList(factors) :=
  57   if zb_operatorp(factors,"^") then
  58     first(factors)
  59   else
  60     if zb_operatorp(factors,"*") then
  61        makelist(factorsList(part(factors,i)),i,1,length(factors))
  62     else
  63        factors;
  64
  65 /* Minimum of a list */
  66 listMin(lst) :=
  67   if length(lst) = 1 then
  68     first(lst)
  69   else
  70     if length(lst) = 2 then
  71        min(first(lst),second(lst))
  72     else
  73        min(first(lst),listMin(rest(lst)));
  74
  75 /* Rising factorial */
  76 rFact(n,k) :=
  77   if k = 0 then
  78     1
  79   else
  80     n*rFact(n+1,k-1);
  81
  82
  83 /* Polynomial coef(deg)*indet^(deg) + coef(deg-1)*indet^(deg) + .... + coef(0)*/
  84 poly(indet,coef,deg) :=
  85   if deg = -1 then
  86     0
  87   else
  88     coef[deg]*indet^deg+poly(indet,coef,deg-1);
  89
  90
  91
  92 /* It creates a polynomial out of a solution of a system of equations */
  93 /* The argument deg must be given the value 0 */
  94 sol2PolyAux(sol,indet,deg) :=
  95   if sol = [] then
  96     0
  97   else
  98     first(sol)*indet^deg + sol2PolyAux(rest(sol),indet,deg+1);
  99
 100 /* Interface macro without auxiliary argument */
 101 sol2Poly(sol,indet) ::=
 102   buildq([sol,indet],
 103          sol2PolyAux(sol,indet,0));
 104
 105 /* It creates a polynomial out of a solution of a system of equations given by equalities */
 106 /* The argument deg must be given the value 0 */
 107 Gsol2PolyAux(sol,indet,deg) :=
 108   if sol = [] then
 109     0
 110   else
 111     rhs(first(sol))*indet^deg + Gsol2PolyAux(rest(sol),indet,deg+1);
 112
 113 /* Interface macro without auxiliary argument */
 114 Gsol2Poly(sol,indet) ::=
 115   buildq([sol,indet],
 116          Gsol2PolyAux(sol,indet,0));
 117
 118
 119 /* Degree of a polynomial (MACRO) */
 120
 121 degree(poly,indet) ::=
 122   buildq([poly,indet],if poly = 0 then -1 else hipow(poly,indet));
 123
 124
 125 /* Leading coefficient of a polynomial */
 126 leadCoeff(poly,indet)::=
 127    buildq([poly,indet],
 128           ratcoeff(poly,indet,hipow(poly,indet)));
 129
 130
 131 /* Leading term of a polynomial */
 132 leadTerm(poly,indet)::=
 133    buildq([poly,indet],
 134           indet^degree(poly,indet));
 135
 136 /* Leading monomial of a polynomial */
 137 leadMono(poly,indet)::=
 138    buildq([poly,indet],
 139           leadCoeff(poly,indet)*leadTerm(poly,indet));
 140
 141 /* Tail of a polynomial */
 142 Tail(poly,indet) ::=
 143    buildq([poly,indet],
 144           poly-leadMono(poly,indet));
 145
 146
 147 /* Given a polynomial, it creates a list of its coefficients */
 148 poly2List(poly,indet) :=
 149   makelist(coeff(poly,indet,i),i,0,degree(poly,indet));
 150
 151
 152 /* Given a list and an array, it creates the corresponding substitution list */
 153 list2subst(list,vector) :=
 154   makelist('vector[i]=list[i],i,0,length(vector));
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165