← Back to branch summary

~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to share/contrib/Zeilberger/shiftQuotient.mac

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2006-10-18 14:52:42 UTC
  • mto: (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20061018145242-vzyrm5hmxr8kiosf
Import upstream version 5.10.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
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(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 op(expr) = binomial then
 
101
      first(expr)!/second(expr)!/(first(expr)-second(expr))!
 
102
    else
 
103
      apply(op(expr),makelist(removeBinomial(part(expr,i)),i,1,length(expr)));
 
104
*/
 
105
 
 
106
removeBinomial(expr) :=
 
107
  block([],
 
108
    
 
109
  if atom(expr) then
 
110
    return(expr)
 
111
  else
 
112
    if op(expr) = binomial then
 
113
      return(first(expr)!/second(expr)!/(first(expr)-second(expr))!)
 
114
    else
 
115
      if op(expr) = "-" then
 
116
        return(-removeBinomial(first(expr)))
 
117
      else
 
118
      
 
119
        return(apply(op(expr),makelist(removeBinomial(part(expr,i)),i,1,length(expr))))
 
120
      
 
121
  );
 
122
 
 
123
shiftFactPoly(expr,k,j) :=
 
124
  if atom(expr) then
 
125
    subst(k+j,k,expr)
 
126
  else
 
127
  if op(expr) = "/" then
 
128
     shiftFactPoly(first(expr),k,j)*shiftFactPoly(second(expr),k,j)^(-1) 
 
129
  else
 
130
  if op(expr) = "*" then
 
131
     product(shiftFactPoly(part(expr,i),k,j),i,1,length(expr)) 
 
132
  else
 
133
  if op(expr) = "^" then
 
134
     if integerp(second(expr)) then 
 
135
         shiftFactPoly(first(expr),k,j)^second(expr)
 
136
     else 
 
137
        0
 
138
  else
 
139
    expand(subst(k+j,k,expr));     
 
140
   
 
141
 
 
142
shiftQuo(expr,k) :=
 
143
  block([sq_res],
 
144
    sq_res : shiftQuoHypCheck(expr,k),
 
145
    if second(sq_res) = HYP then
 
146
      return(first(sq_res))
 
147
    else
 
148
      return(subst(k+1,k,first(sq_res))/first(sq_res))
 
149
    );
 
150
 
 
151
 
 
152
shiftQuoOnlyHyp(expr,k) :=
 
153
  block([sq_res],
 
154
  sq_res : shiftQuoHypCheck(expr,k),
 
155
  if second(sq_res) = HYP then
 
156
    return(first(sq_res))
 
157
  else
 
158
    return([])
 
159
  );
 
160
 
 
161
 
 
162
isPolynomial(expr,k) :=
 
163
  if freeof(k,expr) or expr=k or expr=-k then
 
164
     true
 
165
  else
 
166
    if op(expr) = "^" then
 
167
      freeof(k,second(expr)) and isPolynomial(first(expr),k)
 
168
    else
 
169
      if op(expr) = "-" then
 
170
        isPolynomial(first(expr),k)
 
171
      else
 
172
        if op(expr) = "*" or op(expr)= "+" or op(expr) = "-" then
 
173
           apply("and", makelist(isPolynomial(part(expr,i),k),i,1,length(expr)))
 
174
        else
 
175
           false;
 
176
         
 
177
 
 
178
 
 
179
rationalp(expr,k) :=
 
180
  if not(isPolynomial(expr,k)) and
 
181
     not(op(expr) = "//" and
 
182
                 isPolynomial(expand(num(expr)),k) and
 
183
                 isPolynomial(expand(denom(expr)), k)) then
 
184
     (
 
185
     if op(expr) = "-" then
 
186
       rationalp(first(expr),k)
 
187
     else
 
188
       false
 
189
     )
 
190
  else
 
191
     true;
 
192
 
 
193
shiftQuoHypCheck(expr,k) :=
 
194
  block([xthru_expr,sq_res],
 
195
    xthru_expr : xthru(removeBinomial(expr)),
 
196
    sq_res : shiftQuoHypCheckAux(xthru_expr,k,HYP),
 
197
    
 
198
    if ?ratp(sq_res[1],k) then
 
199
      return(sq_res)
 
200
    else
 
201
      return([sq_res[1],NO_HYP])
 
202
 
 
203
);
 
204
       
 
205
 
 
206
shiftQuoHypCheckAux(expr,k,hyp_flag) :=
 
207
block([sq_num,sq_den,sq_base,sq_exp],
 
208
  if hyp_flag = NO_HYP then
 
209
    return([expr,NO_HYP])
 
210
  else
 
211
  if freeof(k,expr) then
 
212
     return([1,HYP])
 
213
  else
 
214
     if expr = k then
 
215
       return([(k+1)/k,HYP])
 
216
     else
 
217
     if op(expr) = "*" then
 
218
       return(product(shiftQuoHypCheckAux(part(expr,i),k,hyp_flag),
 
219
              i,1,length(expr)))
 
220
     else
 
221
     if op(expr) = "//" then
 
222
       (
 
223
       sq_num : shiftQuoHypCheckAux(first(expr),k,hyp_flag),
 
224
       sq_den : shiftQuoHypCheckAux(second(expr),k,hyp_flag),
 
225
       return([first(sq_num)/first(sq_den),second(sq_num)*second(sq_den)])
 
226
       )
 
227
     else
 
228
     if op(expr) = "^" then
 
229
       (
 
230
      
 
231
       if freeof(k,first(expr)) then
 
232
         sq_base : [first(expr),hyp_flag]
 
233
       else
 
234
         sq_base : shiftQuoHypCheckAux(first(expr),k,hyp_flag),
 
235
       sq_exp : leadCoeff(second(expr),k),
 
236
       if integerp(sq_exp) then
 
237
         return([first(sq_base)^sq_exp,second(sq_base)])
 
238
       else
 
239
         return([first(sq_base)^sq_exp,NO_HYP])
 
240
 
 
241
       )
 
242
     else
 
243
     if op(expr) = "!" then
 
244
        if not(integerp(leadCoeff(first(expr),k))) then
 
245
           return([expr,NO_HYP])
 
246
        else
 
247
          if leadCoeff(first(expr),k)>0 then           
 
248
            return(
 
249
                [product(factor(first(expr)+i),i,1,leadCoeff(first(expr),k)),
 
250
                hyp_flag])
 
251
          else
 
252
             [
 
253
             1/product(factor(first(expr)-i+1), 
 
254
                     i,1,-leadCoeff(first(expr),k)),hyp_flag]
 
255
/*
 
256
     else
 
257
     if op(expr) = binomial then
 
258
        (
 
259
        sq_num : shiftQuoHypCheckAux(factorial(first(expr)),k,hyp_flag),
 
260
        sq_den : shiftQuoHypCheckAux(factorial(first(expr)-second(expr)),
 
261
                                  k,hyp_flag)*   
 
262
                 shiftQuoHypCheckAux(factorial(second(expr)),
 
263
                                  k,hyp_flag),
 
264
        if second(sq_num) = HYP and second(sq_den)=HYP then
 
265
          return([first(sq_num)/first(sq_den),HYP])
 
266
        else
 
267
          return([expr,NO_HYP])
 
268
 
 
269
        )
 
270
*/
 
271
     else
 
272
       if op(expr) = "+" then
 
273
         return([shiftFactPoly(expr,k,1)/expr,hyp_flag])
 
274
       else
 
275
         (
 
276
         print("Unknown operator : ", op(expr)),
 
277
         return([expr,NO_HYP])
 
278
         )
 
279
);
 
280
 
 
281
 
 
282
 
 
283
 
 
284
 
 
285
 
 
286
 
 
287