share/algebra/nusum.mac

   1 /*-*- Mode: macsyma; Package: MAXIMA-*-*/
   2 eval_when([translate,batch,demo,load,loadfile],
   3 dva(var)::=buildq([var],define_variable(var,'var,any)));
   4 dva(%n); dva(%pw); dva(%f); dva(%f1); dva(l%); dva(solvep);
   5 dva(%r); dva(p); dva(%cf);
   6
   7 algebraicp(%1):=catch(
   8            subst("^" = lambda([%1,%2],if not integerp(%2) then throw(true)),%1),
   9            false);
  10 hicoef(x,n):=(x:ratsimp(x,n),coeff(x,n,hipow(x,n)));
  11 genpol(n):=if n < 0 then 0 else concat('%,n)+%n*genpol(n-1);
  12 clist(p):=if 0 = p then []
  13        else cons(ratdisrep(%pw:ratcoef(p,%n,0)),clist(rat((p-%pw)/%n)));
  14 unsum(%g,%n):=if atom(%g) or part(%g,0) # "+"
  15        then factor((num(%g)/subst(%n-1,%n,prodgunch(num(%g),%n,1))
  16                     -denom(%g)/subst(%n-1,%n,prodgunch(denom(%g),%n,1)))
  17                     *(subst(%n-1,%n,num(%g))/denom(%g)))
  18        else map(lambda([x],unsum(x,%n)),%g);
  19 prodflip(%0):=subst(
  20          [nounify('product) = 'product,
  21           'product = lambda([%0,%1,%,%%],1/produ(1/%0,%1,%,%%))],%0);
  22 prodgunch(%0,%n,%2):=
  23   subst([nounify('sin) =
  24     lambda([%0],
  25            sin(subst(%n+%2,%n,%0))
  26           *lambda([trigexpand],
  27                   expand(sin(%0)/sin(subst(%n+%2,%n,%0))))
  28           (true)),
  29          nounify('product) = lambda([%0,%1,%,%3],
  30                                     funmake(nounify('product),
  31                                             [%0,%1,subst(%n+%2,%n,%),
  32                                              subst(%n+%2,%n,%3)])
  33                                    *produ(%0,%1,%,subst(%n+%2,%n,%)-1)
  34                                    /produ(%0,%1,%3+1,subst(%n+%2,%n,%3))),
  35          'binomial = lambda([%0,%1],
  36                             subst(%n+%2,%n,binomial(%0,%1))
  37                            *produ((%1+'%)/(%0+'%),'%,1,subst(%n+%2,%n,%1)-%1)
  38                            *produ((-%1+%0+'%)/(subst(%n+%2,%n,%1)-%1+%0+'%),'%,
  39                                   1,subst(%n+%2,%n,%0-%1)+%1-%0)),
  40          'beta = lambda([%0,%1],
  41                         subst(%n+%2,%n,beta(%0,%1))
  42                        *produ((%0+%1+'%)/(%0+'%),'%,0,ratcoef(%0,%n)*%2-1)
  43                        *produ((%0+%1+%2*ratcoef(%0,%n)+'%)/(%1+'%),'%,0,
  44                               ratcoef(%1,%n)*%2-1)),
  45          "!" = lambda([%0],
  46                       subst(%n+%2,%n,%0)!
  47                         /produ(%0+'%,'%,1,subst(%n+%2,%n,%0)-%0)),
  48          'gamma = lambda([%0],
  49                          gamma(subst(%n+%2,%n,%0))
  50                         /produ(%0+'%-1,'%,1,subst(%n+%2,%n,%0)-%0))],%0);
  51 produ(%0,%1,%,%3):=
  52   lambda([x,y],
  53          if not integerp(y) then funmake(nounify('product),[%0,%1,%,%3])
  54                  else (if y < -1 then 1/produ(%0,%1,%3+1,%-1)
  55                            else (for i from 0 thru y do x:x*subst(i+%,%1,%0),x)))(
  56       1,ratsimp(%3-%));
  57 nusum(%a,%n,%l,%h):=block([mapprint:false,programmode:true,solvenullwarn:false],
  58       nusuml(%a,%n,%l,%h,[])[1]);
  59 nusum(%a,%n,%l,%h):=block([mapprint:false,programmode:true,solvenullwarn:false],
  60       first(nusuml(%a,%n,%l,%h,[])));
  61 funcsolve(%a,%f):=block([mapprint:false,programmode:true,solvenullwarn:false],
  62           funcsol(%a,%f,[])[1]);
  63 dimsum(%cl):=block([ratfac:false,%cd,%pt,%pw],
  64        %cd:map(lambda([x],hipow(ratsimp(x)+1/%n,%n)),
  65                [%cl[1]+%cl[2],%cl[1]-%cl[2],%cl[3]]),%cd[1]:max(%cd[1],%cd[2]-1),
  66        inpart(subst(solvep:solve(clist(subst(
  67                                         %pt
  68                                          :funmake('lambda,
  69                                                   [[%n],
  70                                                    genpol(
  71                                                     if
  72                                                      %cd[1] < %cd[2]
  73                                                        and integerp(
  74                                                        %pw
  75                                                         :subst(
  76                                                          solve(
  77                                                           ratcoef(
  78                                                            %cl[1]*(%n+'%)
  79                                                             +%cl[2]*%n,%n,%cd[2]),
  80                                                           '%),'%))
  81                                                       then max(%pw,%cd[3]-%cd[1])
  82                                                       else %cd[3]-%cd[1])]),
  83                                         inpart(%f,0),%cl . [%f1,%f,1])),
  84                                  append(if listp(l%) then l% else l%:[l%],
  85                                         clist(inpart(%pt,2)))),
  86                     (l%:map("=",l%,subst(solvep,l%)),%pt)),2));
  87
  88 ratsolve(p,x):=apply('append,
  89          maplist(lambda([p],
  90                         if hipow(p,x) = 1 or subst(0,x,p) = 0
  91                             then solve(subst(lambda([x,y],x),"^",p),x) else []),
  92                  2*factor(p)^2));
  93 prodshift(%0,%2):=subst(
  94           [nounify('product) = 'product,
  95            'product = lambda([%0,%1,%,%3],produ(subst(%1-%2,%1,%0),%1,%+%2,%3+%2))],
  96           %0);
  97
  98 rforn(%3):=block([y:gcd(%r[2],subst(%n-%3,%n,%r[1])),ratfac:true],
  99       p:p*produ(y,%n,%n,%n+%3-1),%r:ratsimp(%r/[subst(%n+%3,%n,y),y]))$
 100
 101 /*  Note that this business of w/[a,b] is very sensitive to evaluation.  It
 102   is assumed that w is evaluated first but sometimes this doesn't happen
 103   resulting in wrong answer eg w/[a,b]==> [w/a,w/b] and then eval w:[1,2] bad!!*/
 104 /*rforn(%3):=block([y:gcd(%r[2],subst(%n-%3,%n,%r[1])),ratfac:true,vv],
 105       p:p*produ(y,%n,%n,%n+%3-1),%r[1]:ratsimp(%r[1]/subst(%n+%3,%n,y)),%r[2]:ratsimp(%r[2]/y));*/
 106
 107 rform(%r):=if ?ratp(%r[1]/%r[2],%n)
 108        then (if algebraicp(%r) then (gcd:'red,algebraic:true),
 109              block([p:1],rforn(1),
 110                    for %3 in ratsolve(resultant(%r[1],subst(%n+'%,%n,%r[2]),%n),'%)
 111                        do
 112                        (if integerp(%3:subst([%3],'%)) and %3 > 0 then rforn(%3)),
 113                    [p,%r[1]/%r[2]]))
 114        else error(%r[1]/%r[2],"non-rational term ratio to nusum");
 115 nusuml(%a,%n,%l,%h,l%):=
 116   if %a = 0
 117     then [0]
 118     else
 119       block([solvep:false,modulus:false,rv:ratvars,prodhack:true,
 120              ratfac:true,gcd:'spmod,algebraic:false,ratalgdenom:true,
 121              matrix_element_mult:"*",
 122              dispflag:false,maperror:false,%f:funmake('%f,[%n]),
 123              %f1:funmake('%f,[%n+1])],
 124         ratvars(%n),
 125         if [] # errcatch
 126           (%cf:dimsum
 127             ([-num(
 128                   (%r
 129                     :rform(
 130                           lambda([%a],[num(%a),denom(%a)])
 131                           (factor(
 132                                  subst(%n+1,%n,%a)
 133                                  /prodgunch(%a,%n,1)))))[
 134                                    2]),subst(%n-1,%n,denom(%r[2])),%r[1]]))
 135            and [] # solvep
 136              then cons((%f:prodgunch(num(%a),%n,1)/denom(%a),
 137                         %f1:ratsimp(radcan(%cf)),apply('ratvars,rv),
 138                         %f1:subst(lambda([%0,%1,%,%3],produ(%0,%1,%,%3)),
 139                                   nounify('product),
 140                                   factor(
 141                                         subst(%h,%n,
 142                                               factor(
 143                                                     num(%r[2])*%f
 144                                                     *subst(%n+1,%n,%f1)
 145                                                     /%r[1])))
 146                                  -factor(
 147                                         subst(%l,%n,
 148                                               factor(
 149                                                     %a*%f1
 150                                                     *subst(%n-1,%n,denom(%r[2]))
 151                                                     /%r[1])))),
 152                         if ?ratp(%a,%n) then factor(%f1) else %f1),l%)
 153                                         else (apply('ratvars,rv),[apply('sum,[%a,%n,%l,%h])]));
 154
 155 funcsol(%a,%f,l%):=block(
 156   [ratfac:true,maperror:false,linenum:linenum,dispflag:false,%f1,%cl,%cm,%n:inpart(%f,1)],
 157         %f1:subst(%n = %n+1,%f),
 158         %cl:factor(augcoefmatrix([%a:num(xthru(lhs(%a)-rhs(%a)))],[%f1,%f])[1]),
 159         %cm:rform(rest(%cl,-1)),%cm[2]:ratsimp(subst(%n+1,%n,%cm[1])/%cm[1]),
 160         append(errcatch(%f = factor(dimsum(
 161                                      ratsimp(
 162                                       [%cl[1]/num(%cm[2]),%cl[2]/denom(%cm[2]),
 163                                        %cm[1]*%cl[3]/denom(%cm[2])]))
 164                                      /%cm[1])),l%));
 165
 166