module sum2; % Auxiliary package for summation in finite terms.

% Authors: K.Yamamoto, K.Kishimoto & K.Onaga  Hiroshima Univ.
% Modified by: F.Kako                         Hiroshima Univ.
%         Fri Sep. 19, 1986
%         Mon Sep. 7, 1987   added PROD operator (by F. Kako)
% e-mail: kako@kako.math.sci.hiroshima-u.ac.jp
%         or D52789%JPNKUDPC.BITNET

% Usage:
%         sum(expression,variable[,lower[,upper]]);
%         lower and upper are optionals.

%         prod(expression,variable[,lower[,upper]]);
%         returns product of expression.

fluid '(!*trsum);               % trace switch;

fluid '(sum_last_attempt_rules!*);

switch trsum;

symbolic procedure simp!-sum0(u,y);
  begin scalar v,upper,lower,lower1,dif;
      if not atom cdr y then <<
          lower := cadr y;
          lower1 := if numberp lower then lower - 1
                     else list('plus,lower,-1);
          upper := if not atom cddr y then caddr y else car y;
          dif := addsq(simp!* upper, negsq simp!* lower);
          if denr dif = 1 then
              if null numr dif
                then return subsq(u,list(!*a2k car y . upper))
                else if fixp numr dif then dif := numr dif
                else dif := nil
           else dif := nil;
          if dif and dif <= 0 then return nil ./ 1;
          if atom cddr y then upper := nil>>;
      v := !*a2k car y;
      return simp!-sum1(u,v,y,upper,lower,lower1,dif)
   end;

symbolic procedure simp!-sum1(u,v,y,upper,lower,lower1,dif);
   begin scalar w,lst,x,z,flg;
      lst := sum!-split!-log(u,v);
      w := car lst;
      lst := cdr lst;
      u := nil ./ 1;
    a:
      if null w then go to b;
      x := multsq(caar w,
                  simp!-prod1(cdar w,v,y,upper,lower,lower1,dif));
      u := addsq(u,simp!* list('log, prepsq x));
      w := cdr w;
      go to a;
    b:
       if null lst then return u;
            flg := nil;
            z := car lst;
            if !*trsum then <<
               prin2!* "Summation ";sqprint z;prin2!*  " w.r.t ";
               xprinf(!*k2f v,nil,nil);terpri!* t >>;
%           z := reorder numr z ./ reorder denr z;
            w := sum!-sq(z,v);
            if w = 'failed then <<
                if !*trsum then <<
                     prin2!* "UMM-SQ failed. Trying SUM-TRIG";
                     terpri!* t>>;
                w := sum!-trig(z,v);
                if w = 'failed then <<
                     if !*trsum then <<
                         prin2!* "SUM-TRIG failed.";
                         terpri!* t>>;
                     w := sum!-unknown(z,v,y,lower,dif);
                     flg := car w;
                     w := cdr w>>>>;
            if !*trsum then <<
                 prin2!* "Result  = "; sqprint w; terpri!* t >>;
            if flg then goto c;
            if upper then
                w := addsq(subsq(w,list(v . upper)),
                           negsq subsq(w,list(v . lower1)))
             else if lower then
                w := addsq(w , negsq subsq(w, list(v . lower1)));
   c:
            u := addsq(u,w);
            lst := cdr lst;
            goto b
   end;


%*********************************************************************
%       Case of trigonometric or other functions
%       Trigonometric functions are expressed in terms of exponetials.
%       Pattern matching to get the summation in closed form.
%********************************************************************;

global '(!*trig!-to!-exp);              % variable to indicate
                                        % that the expression contains
                                        % some trig. functions.

symbolic procedure sum!-trig(u,v);
   begin scalar lst,w;  % z;
      !*trig!-to!-exp := nil;           % trig. to exponential.
      u := trig!-to!-expsq(u,v);
      if not !*trig!-to!-exp then return 'failed;
      lst := sum!-term!-split(u,v);
      u := nil ./ 1;
   a:
      if null lst then return exp!-to!-trigsq u;
%     z := reorder numr car lst ./ reorder denr car lst;
%     w := sum!-sq(z,v);
      w := sum!-sq(car lst,v);
      if w = 'failed then return 'failed;
%      w := exp!-to!-trigsq w;     % exponential to trig. function.
      u := addsq(u,w);
      lst := cdr lst;
      goto a
   end;

   sum_last_attempt_rules!* :=
    algebraic  <<
    { sum(~f + ~g,~n,~anf,~ende) => sum(f,n,anf,ende) +
                                       sum(g,n,anf,ende)
                when or (part(sum(f,n,anf,ende),0) neq sum ,
                         part(sum(g,n,anf,ende),0) neq sum ),
      sum((~f+~g)/~dd,~n,~anf,~ende) => sum(f/dd,n,anf,ende) +
                                        sum(g/dd,n,anf,ende)
                when or (part(sum(f/dd,n,anf,ende),0) neq sum ,
                         part(sum(g/dd,n,anf,ende),0) neq sum ),
      sum(~c*~f,~n,~anf,~ende) => c* sum(f,n,anf,ende)
                         when freeof(c,n) and c neq 1,
      sum(~c/~f,~n,~anf,~ende) => c* sum(1/f,n,anf,ende)
                         when freeof(c,n) and c neq 1,
      sum(~c*~f/~g,~n,~anf,~ende) => c* sum(f/g,n,anf,ende)
                          when freeof(c,n) and c neq 1} >>$


symbolic procedure sum!-unknown(u,v,y,lower,dif);
   begin scalar z,w;
     if null dif then <<
          z := 'sum . (prepsq u . list car y);
          if w := opmtch z then return (nil . simp w)
           else if null cdr y then return (t . mksq(z,1));
          z := 'sum . (prepsq u . y);
          let sum_last_attempt_rules!*;
          w:= opmtch z;
          rule!-list (list sum_last_attempt_rules!*,nil);
          return (t . if w then simp w else mksq(z,1))>>;
%         return (t . if w := opmtch z then simp w else mksq(z,1))>>;
     z := nil ./ 1;
  a:
     if dif < 0 then return (t . z);
     z := addsq(z,subsq(u,list(v . list('plus,lower,dif))));
     dif := dif - 1;
     goto a
   end;


%*********************************************************************
%       Summation by Gosper's algorithm.
%********************************************************************;

symbolic procedure sum!-sq(u,v);
   %Argument U : expression of s-q;
   %         V : kernel;
   %value        : expression of sq (result of summation.);
   begin scalar gn,fn,pn,rn,qn,z,k,x;
      if null numr u then return nil ./ 1;
      x := setkorder list v;
      z := reorder numr u;
      u := z ./ reorder denr u;
      if !*trsum then <<
            prin2t " *** Summation by Gosper's algorithm ***";
            prin2!* "    A(n) = "; sqprint u;terpri!* t;
            terpri!* t>>;
      if domainp z or not (mvar z eq v) or red z then
            <<pn := 1; z := u>>
       else <<pn := !*p2f lpow z;
              z := lc z ./ denr u>>;
      z := quotsq(z,nsubsq(z,v, - 1));
      gn := gcdf!*(numr z,denr z);
      if !*trsum then <<
         prin2!* "A(n)/A(n-1) = ";sqprint z;terpri!* t;
         prin2!* "GN = ";xprinf(gn,nil,nil);terpri!* t>>;
      qn := quotf!*(numr z, gn);
      rn := quotf!*(denr z, gn);
      if nonpolyp(qn,v) or nonpolyp(rn,v) then go to fail;
      if !*trsum then <<
         prin2!* "Initial qn, rn and pn are "; terpri!* t;
         prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t;
         prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t;
         prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>;
      k := compress explode '!+j;
      z := integer!-root(resultant(qn,nsubsf(rn,v,k),v),k);
      if !*trsum then <<
         prin2 "Root of resultant(q(n),r(n+j)) are "; prin2t z >>;
      while z do <<
            k := car z;
            gn := gcdf!*(qn,nsubsf(rn,v,k));
            qn := quotf!*(qn,gn);
            rn := quotf!*(rn,nsubsf(gn,v, -k));
            while (k := k - 1)>=0 do
               pn := multf(pn,nsubsf(gn,v, -k));
            z := cdr z>>;
      if !*trsum then <<
         prin2!* "Shift free  qn, rn and pn are";terpri!* t;
         prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t;
         prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t;
         prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>;
      qn := nsubsf(qn,v,1);
      if (k := degree!-bound(pn,addf(qn,rn),addf(qn,negf rn),v)) < 0
         then go to fail;
      if !*trsum then <<
         prin2 "DEGREE BOUND is "; prin2t k >>;
      if not(fn := solve!-fn(k,pn,qn,rn,v)) then go to fail;
      if !*trsum then <<
         prin2!* "FN = ";sqprint fn;terpri!* t >>;
      u := multsq(multsq(qn ./ 1,fn), multsq(u, 1 ./ pn));
      z := gcdf!*(numr u, denr u);
      u := quotf!*(numr u, z) ./ quotf!*(denr u,z);
      if !*trsum then <<
            prin2t " *** Gosper's algorithm completed ***";
            prin2!* "    S(n) = "; sqprint u;terpri!* t;
            terpri!* t>>;
      setkorder x;
      return (reorder numr u ./ reorder denr u);
  fail:
      if !*trsum then <<
            prin2t " *** Gosper's algorithm failed ***";
            terpri!* t>>;
      setkorder x;
      return 'failed
   end;

%*********************************************************************
%       integer root isolation
%********************************************************************;

symbolic procedure integer!-root(u,v);
% Produce a list of all positive integer root of U;
   begin scalar x,root,n,w;
        x := setkorder list v;
        u := reorder u;
        if domainp u or not(mvar u eq v) then go to a;
        u := numr cancel(u ./ lc u);
        w := u;                                 % get trailing term;
        while not domainp w and mvar w eq v and cdr w do
            w := cdr w;
        if (n := degr(w,v)) > 0 then <<
            w := lc w;
            while n > 0 do <<
               root := 0 . root;
               n := n - 1>>>>;
        n := dfactors lowcoef w;                % factor tail coeff.
        w := (v . 1) . 1;
        while n do <<
            if not testdivide(u,v,car n) then <<
                root := car n . root;
                u := quotf!*(u, (w . - car n))>>
             else n := cdr n>>;
  a:
        setkorder x;
        return root
   end;


symbolic procedure lowcoef u;
   begin scalar lst,m;
        lst := dcoefl u;
        m := 0;
  a:
        if null lst then return m;
        m := gcdn(m,car lst);
        if m = 1 then return 1;
        lst := cdr lst;
        go to a
   end;


symbolic procedure dcoefl u;
   if domainp u then if fixp u then list abs u else nil
    else nconc(dcoefl lc u , dcoefl red u);

symbolic procedure testdivide(u,v,n);
% Evaluate U at integer point (V = N);
   begin scalar x;
 a:
       if domainp u or not(mvar u eq v) then return addf(u,x);
       x := addf(multd(expt(n,ldeg u),lc u),x);
       if (u := red u) then go to a;
       return x
   end;


%*********************************************************************
%********************************************************************;

symbolic procedure degree!-bound(pn,u,v,kern);
% degree bound for fn;
%       u: q(n+1) + r(n);
%       v: q(n+1) - r(n);
   begin scalar lp,l!+, l!-, x,m,k;
      x := setkorder list kern;
      u := reorder u;
      v := reorder v;
      pn := reorder pn;
      l!+ := if u then degr(u,kern) else -1;
      l!- := if v then degr(v,kern) else -1;
      lp := if pn then degr(pn,kern) else -1;
      if l!+ <= l!- then <<k := lp - l!-;go to  a>>;
      k := lp - l!+ + 1;
      if l!+ > 0 then u := lc u;
      if l!- > 0 then v := lc v;
      if l!+ = l!- + 1 and fixp(m := quotf1(multd(-2,v),u)) then
        k := max(m,k)
       else if lp = l!- then k := max(k,0);
   a:
      setkorder x;
      return k
   end;

%*********************************************************************
%       calculate polynomial f(n) such that
%       p(n) - q(n+1)*f(n) + r(n)*f(n-1) = 0;
%********************************************************************;

symbolic procedure solve!-fn(k,pn,qn,rn,v);
   begin scalar i,fn,x,y,z,u,w,c,clst,flst;
      c := makevar('c,0);
      clst := list c;
      fn := !*k2f c;
      i := 0;
      while (i := i + 1) <= k do <<
         c := makevar('c,i);
         clst := c . clst;
         fn := ((v . i) . !*k2f c) . fn>>;
      z :=
       addf(pn,
            addf(negf multf(qn,fn),
                 multf(rn,nsubsf(fn,v, - 1))));
      x := setkorder (v . clst);
      z := reorder z;
      c := clst;
      if !*trsum then <<
         prin2!* "C Equation is";terpri!* t;
         xprinf(z,nil,nil);terpri!* t >>;
  a:
      if domainp z or
         domainp (y := if mvar z eq v then lc z else z) then
           go to fail;
      w := mvar y;
      if not(w memq clst) then go to fail;
      if !*trsum then <<
         prin2!* "C Equation to solve is ";xprinf(y,nil,nil);terpri!* t;
         prin2!* " w.r.t ";xprinf(!*k2f w,nil,nil);terpri!* t >>;
      u := gcdf!*(red y , lc y);
      u := quotf!*(negf red y, u) ./ quotf!*(lc y, u);
      flst := (w . u) . flst;
      z := subst!-cn(z,w,u);
      if !*trsum then <<
         xprinf(!*k2f w,nil,nil);prin2!* " := ";sqprint u;terpri!* t >>;
      clst := deleteq(clst,w);
      if z then go to a;
      setkorder c;
      fn := reorder fn;
      u := 1;
      while not domainp fn and mvar fn memq c do <<
         w := mvar fn;
         z := atsoc(w,flst);
         fn := subst!-cn(fn,w,if z then cdr z);
         if z then u := multf(u,denr cdr z);
         z := gcdf!*(fn,u);
         fn := quotf!*(fn,z);
         u := quotf!*(u,z)>>;
      setkorder x;
      return cancel(reorder fn ./ reorder u);
 fail:
     if !*trsum then <<
        prin2t "Fail to solve C equation.";
        prin2!* "Z := ";xprinf(z,nil,nil);terpri!* t >>;
     setkorder x;
     return nil
   end;


symbolic procedure subst!-cn(u,v,x);
   begin scalar z;
      z := setkorder list v;
      u := reorder u;
      if not domainp u and mvar u eq v then
         if x then u := addf(multf(lc u,reorder numr x),
                             multf(red u,reorder denr x))
          else u := red u;
      setkorder z;
      return reorder u
   end;



symbolic procedure makevar(id,n);
  compress nconc(explode id, explode n);

symbolic procedure deleteq(u,x);
   if null u then nil
    else if car u eq x then cdr u
    else car u . deleteq(cdr u, x);


symbolic procedure nsubsf(u,kern,i);
% ARGUMENT U : expression of sf;
%          KERN : kernel;
%          I : integer or name of integer variable;
% value        : expression of sf;
   begin scalar x,y,z,n;
      if null i or i = 0 then return u;
      x := setkorder list kern;
      u := reorder u;
      y := addf(!*k2f kern,
                if fixp i then i else !*k2f i);
      z := nil;
   a:
      if domainp u or not(mvar u eq kern) then goto b;
      z := addf(z,lc u);
      n := degr(u,kern) - degr(red u,kern);
      u := red u;
   a1:
      if n <= 0 then goto a;
      z := multf(z,y);
      n := n - 1;
      go to a1;
   b:
      z := addf(z,u);
      setkorder x;
      return reorder z
   end;


symbolic procedure nsubsq(u,kern,i);
% ARGUMENT U : expression of sq;
%       KERN : kernel;
%          I : integer or name of integer variable;
% value        : expression of sq;
   subsq(u,list(kern . list('plus, kern, i)));


%*********************************************************************
%       dependency check
%********************************************************************;

symbolic procedure nonpolyp(u,v);
% check U is not a polynomial of V;
   if domainp u then nil
    else (not(mvar u eq v) and depend!-p(mvar u,v))
         or nonpolyp(lc u,v) or nonpolyp(red u,v);


symbolic procedure depend!-sq(u,v);
   depend!-f(numr u,v) or depend!-f(denr u,v);


symbolic procedure depend!-f(u,v);
   if domainp u then nil
    else depend!-p(mvar u,v) or
         depend!-f(lc u,v) or
         depend!-f(red u,v);


symbolic procedure depend!-p(u,v);
   if u eq v then t
    else if atom u then nil
    else if not atom car u then depend!-f(u,v)
    else if car u eq '!*sq then depend!-sq(cadr u, v)
    else depend!-l(cdr u, v);

symbolic procedure depend!-l(u,v);
   if null u then nil
    else if depend!-sq(simp car u, v) then t
    else depend!-l(cdr u,v);


%*********************************************************************
%       term splitting
%********************************************************************;

symbolic procedure sum!-term!-split(u,v);
   begin scalar y,z,klst,lst,x;
      x := setkorder list v;
      z := qremf(reorder numr u, y := reorder denr u);
      klst := kern!-list(car z,v);
      lst := termlst(car z, 1 ./ 1, klst);
      klst := kern!-list(cdr z,v);
      if depend!-f(y,v) then klst := deleteq(klst,v);
      lst := append(lst, termlst(cdr z, 1 ./ y, klst));
      setkorder x;
      return lst
   end;


symbolic procedure kern!-list(u,v);
% Returns list of kernels that depend on V;
   begin scalar x;
      for each j in kernels u do if depend!-p(j,v) then x := j . x;
      return x
   end;

symbolic procedure termlst(u,v,klst);
   begin scalar x,kern,lst;
      if null u then return nil
       else if null klst or domainp u then
                 return list multsq(!*f2q u,v);
      kern := car klst;
      klst := cdr klst;
      x := setkorder list kern;
      u := reorder u;
      v := reorder(numr v) ./ reorder(denr v);
      while not domainp u and mvar u eq kern do <<
         lst := nconc(termlst(lc u, multsq(!*p2q lpow u, v),klst),lst);
         u := red u>>;
      if u then lst := nconc(termlst(u,v,klst),lst);
      setkorder x;
      return lst
   end;


%*********************************************************************
%       Express trigonometric functions (such as sin, cos ..)
%       by exponentials.
%********************************************************************;

symbolic procedure trig!-to!-expsq(u,v);
   multsq(trig!-to!-expf(numr u,v),
          invsq trig!-to!-expf(denr u,v));


symbolic procedure trig!-to!-expf(u,v);
   if domainp u then u ./ 1
    else addsq(multsq(trig!-to!-expp(lpow u,v),
                      trig!-to!-expf(lc u,v)),
               trig!-to!-expf(red u,v));


symbolic procedure trig!-to!-expp(u,v);
   begin scalar !*combineexpt,w,x,z,n,wi;
      % We don't want to combine expt terms here, since the code
      % depends on the terms being separate.
      n := cdr u;          % integer power;
      z := car u;          % main variable;
      if atom z or not atom (x := car z) or not depend!-p(z,v) then
           return !*p2q u;
      if x memq '(sin cos tan sec cosec cot) then <<
          !*trig!-to!-exp := t;
          w := multsq(!*k2q 'i, simp!* cadr z);
          w := simp!* list('expt,'e, mk!*sq w);
%         W := SIMP LIST('EXPT,'E, 'TIMES . ( 'I . CDR Z));
          wi := invsq w;
          if x eq 'sin then
              w := multsq(addsq(w ,negsq wi),
                          1 ./ list(('i .** 1) .* 2))
           else if x eq 'cos then
              w := multsq(addsq(w, wi), 1 ./ 2)
           else if x eq 'tan then
              w := multsq(addsq(w,negsq wi),
                          invsq addsq(w,wi))
           else if x eq 'sec then
              w := multsq(2 ./ 1, invsq addsq(w, wi))
           else if x eq 'cosec then
              w := multsq(list(('i .** 1) .* 2),
                          invsq addsq(w, negsq wi))
           else
              w := multsq(addsq(w, wi),
                          invsq addsq(w, negsq wi))
                     >>
       else if x memq '(sinh cosh tanh sech cosech coth) then <<
          !*trig!-to!-exp := t;
          w := simp!* list('expt,'e,cadr z);
          wi := invsq w;
          if x eq 'sinh then
              w := multsq(addsq(w,negsq wi), 1 ./ 2)
           else if x eq 'cosh then
              w := multsq(addsq(w,wi), 1 ./ 2)
           else if x eq 'tanh then
              w := multsq(addsq(w,negsq wi),
                          invsq addsq(w,wi))
           else if x eq 'sech then
              w := multsq(2 ./ 1, invsq addsq(w, wi))
           else if x eq 'cosech then
              w := multsq(2 ./ 1, invsq addsq(w, negsq wi))
           else
              w := multsq(addsq(w,wi),
                          invsq addsq(w, negsq wi))
                 >>
       else return !*p2q u;
      return exptsq(w,n)
   end;

%*********************************************************************
%       Inverse of trig!-to!-exp.
%       Express exponentials in terms of trigonometric functions
%                (sin, cos, sinh and cosh)
%                               Wed Dec. 17, 1986 by F. Kako;
%********************************************************************;

symbolic procedure exp!-to!-trigsq u;
   multsq(exp!-to!-trigf numr u,
          invsq exp!-to!-trigf denr u);

symbolic procedure exp!-to!-trigf u;
   begin scalar v,v1,x,y,n;
      u := termlst1(u,1,nil ./1);
      v := nil;
   a:
      if null u then go to b;
      x := caar u;
      y := cdar u;
      u := cdr u;
   a1:
      if u and y = cdar u then <<
         x := addf(x,caar u);
         u := cdr u;
         go to a1>>;
      v := (x . y) . v;
      go to a;
   b:
      v1 := reverse v;
      n := length v;
      u := nil ./ 1;
   c:
      if n = 0 then return u
       else if n = 1 then
          return addsq(u,
                       multsq(!*f2q caar v,
                              simp!* list('expt,'e,mk!*sq cdar v)));
      u := addsq(u,exp!-to!-trigl(caar v1,caar v,cdar v1,cdar v));
      v := cdr v;
      v1 := cdr v1;
      n := n - 2;
      go to c
   end;

symbolic procedure exp!-to!-trigl(a,b,c,d);
%       A*E**C + B*E**D
%               -->
%       ((A+B)*COSH((C-D)/2)+(A-B)*SINH((C-D)/2))*E**((C+D)/2);
%       A, B: sf;
%       C, D: sq;
   begin scalar x,y,z;
      x := !*f2q addf(a,b);
      y := !*f2q addf(a, negf b);
      z := multsq(addsq(c,negsq d), 1 ./ 2);
      z := real!-imag!-sincos z;
      return multsq(simp!* list('expt,'e,
                                mk!*sq multsq(addsq(c,d),1 ./ 2)),
                    addsq(multsq(x, cdr z),
                          multsq(y, car z)))
   end;


symbolic procedure termlst1(u,v,w);
%ARGUMENT U : sf;
%         V : sf;
%         W : sq;
%value      : list of (sf . sq);
   begin scalar x,y;
      if null u then return nil
       else if domainp u then return list (multf(u,v) . w);
      x := mvar u;
      y := if atom x or not(car x eq 'expt) or not(cadr x eq 'e) then
              termlst1(lc u,multf(!*p2f lpow u,v),w)
            else termlst1(lc u,v,addsq(w,multsq(simp!* caddr x,
                                                ldeg u ./ 1)));
      return nconc(y,termlst1(red u,v,w))
   end;


% These can be found in Abramowitz-Stegun (27.8.6 Summable Series), and
% were suggested by Winfried Neun.

algebraic;

let {sum(sin(~n*~tt)/n,~n,1,infinity) => (pi - tt)/2,
     sum(sin(~n*~tt)/(~n)^3,~n,1,infinity) =>
             pi^2*tt/6 - pi*tt^2/4 + tt^3/12,
     sum(sin(~n*~tt)/(~n)^5,~n,1,infinity) =>
             pi^4*tt/90 - pi^2*tt^3/36 + pi*tt^4/48-tt^5/240}$

let {sum(cos(~n*~tt)/(~n),~n,1,infinity) =>  -log(2*sin(tt/2)),
     sum(cos(~n*~tt)/(~n)^2,~n,1,infinity) =>
             pi^2/6 - pi*tt/2 + tt^2/4,
     sum(cos(~n*~tt)/(~n)^4,~n,1,infinity) =>
             pi^4/90 - pi^2*tt^2/12 + pi*tt^3/12-tt^4/48}$

endmodule;

end;
