1module sum2; % Auxiliary package for summation in finite terms.
2
3% Authors: K.Yamamoto, K.Kishimoto & K.Onaga  Hiroshima Univ.
4% Modified by: F.Kako                         Hiroshima Univ.
5%         Fri Sep. 19, 1986
6%         Mon Sep. 7, 1987   added PROD operator (by F. Kako)
7% e-mail: kako@kako.math.sci.hiroshima-u.ac.jp
8%         or D52789%JPNKUDPC.BITNET
9
10% Redistribution and use in source and binary forms, with or without
11% modification, are permitted provided that the following conditions are met:
12%
13%    * Redistributions of source code must retain the relevant copyright
14%      notice, this list of conditions and the following disclaimer.
15%    * Redistributions in binary form must reproduce the above copyright
16%      notice, this list of conditions and the following disclaimer in the
17%      documentation and/or other materials provided with the distribution.
18%
19% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
21% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
23% CONTRIBUTORS
24% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30% POSSIBILITY OF SUCH DAMAGE.
31%
32
33
34% Usage:
35%         sum(expression,variable[,lower[,upper]]);
36%         lower and upper are optionals.
37
38%         prod(expression,variable[,lower[,upper]]);
39%         returns product of expression.
40
41fluid '(!*trsum);               % trace switch;
42
43fluid '(sum_last_attempt_rules!*);
44
45switch trsum;
46
47symbolic procedure simp!-sum0(u,y);
48  begin scalar v,upper,lower,lower1,dif;
49      if not atom cdr y then <<
50          lower := cadr y;
51          lower1 := if numberp lower then lower - 1
52                     else list('plus,lower,-1);
53          upper := if not atom cddr y then caddr y else car y;
54          dif := addsq(simp!* upper, negsq simp!* lower);
55          if denr dif = 1 then
56              if null numr dif
57                then return subsq(u,list(!*a2k car y . upper))
58                else if fixp numr dif then dif := numr dif
59                else dif := nil
60           else dif := nil;
61          if dif and dif <= 0 then return nil ./ 1;
62          if atom cddr y then upper := nil>>;
63      v := !*a2k car y;
64      return simp!-sum1(u,v,y,upper,lower,lower1,dif)
65   end;
66
67symbolic procedure simp!-sum1(u,v,y,upper,lower,lower1,dif);
68   begin scalar w,lst,x,z,flg;
69      lst := sum!-split!-log(u,v);
70      w := car lst;
71      lst := cdr lst;
72      u := nil ./ 1;
73    a:
74      if null w then go to b;
75      x := multsq(caar w,
76                  simp!-prod1(cdar w,v,y,upper,lower,lower1,dif));
77      u := addsq(u,simp!* list('log, prepsq x));
78      w := cdr w;
79      go to a;
80    b:
81       if null lst then return u;
82            flg := nil;
83            z := car lst;
84            if !*trsum then <<
85               prin2!* "Summation ";sqprint z;prin2!*  " w.r.t ";
86               xprinf(!*k2f v,nil,nil);terpri!* t >>;
87%           z := reorder numr z ./ reorder denr z;
88            w := sum!-sq(z,v);
89            if w = 'failed then <<
90                if !*trsum then <<
91                     prin2!* "UMM-SQ failed. Trying SUM-TRIG";
92                     terpri!* t>>;
93                w := sum!-trig(z,v);
94                if w = 'failed then <<
95                     if !*trsum then <<
96                         prin2!* "SUM-TRIG failed.";
97                         terpri!* t>>;
98                     w := sum!-unknown(z,v,y,lower,dif);
99                     flg := car w;
100                     w := cdr w>>>>;
101            if !*trsum then <<
102                 prin2!* "Result  = "; sqprint w; terpri!* t >>;
103            if flg then goto c;
104            if upper then
105                w := addsq(subsq(w,list(v . upper)),
106                           negsq subsq(w,list(v . lower1)))
107             else if lower then
108                w := addsq(w , negsq subsq(w, list(v . lower1)));
109   c:
110            u := addsq(u,w);
111            lst := cdr lst;
112            goto b
113   end;
114
115
116%*********************************************************************
117%       Case of trigonometric or other functions
118%       Trigonometric functions are expressed in terms of exponetials.
119%       Pattern matching to get the summation in closed form.
120%********************************************************************;
121
122global '(!*trig!-to!-exp);              % variable to indicate
123                                        % that the expression contains
124                                        % some trig. functions.
125
126symbolic procedure sum!-trig(u,v);
127   begin scalar lst,w;  % z;
128      !*trig!-to!-exp := nil;           % trig. to exponential.
129      u := trig!-to!-expsq(u,v);
130      if not !*trig!-to!-exp then return 'failed;
131      lst := sum!-term!-split(u,v);
132      u := nil ./ 1;
133   a:
134      if null lst then return exp!-to!-trigsq u;
135%     z := reorder numr car lst ./ reorder denr car lst;
136%     w := sum!-sq(z,v);
137      w := sum!-sq(car lst,v);
138      if w = 'failed then return 'failed;
139%      w := exp!-to!-trigsq w;     % exponential to trig. function.
140      u := addsq(u,w);
141      lst := cdr lst;
142      goto a
143   end;
144
145   sum_last_attempt_rules!* :=
146    algebraic  <<
147    { sum(~f + ~g,~n,~anf,~ende) => sum(f,n,anf,ende) +
148                                       sum(g,n,anf,ende)
149                when or (part(sum(f,n,anf,ende),0) neq sum ,
150                         part(sum(g,n,anf,ende),0) neq sum ),
151      sum((~f+~g)/~dd,~n,~anf,~ende) => sum(f/dd,n,anf,ende) +
152                                        sum(g/dd,n,anf,ende)
153                when or (part(sum(f/dd,n,anf,ende),0) neq sum ,
154                         part(sum(g/dd,n,anf,ende),0) neq sum ),
155      sum(~c*~f,~n,~anf,~ende) => c* sum(f,n,anf,ende)
156                         when freeof(c,n) and c neq 1,
157      sum(~c/~f,~n,~anf,~ende) => c* sum(1/f,n,anf,ende)
158                         when freeof(c,n) and c neq 1,
159      sum(~c*~f/~g,~n,~anf,~ende) => c* sum(f/g,n,anf,ende)
160                          when freeof(c,n) and c neq 1} >>$
161
162
163symbolic procedure sum!-unknown(u,v,y,lower,dif);
164   begin scalar z,w;
165     if null dif then <<
166          z := 'sum . (prepsq u . list car y);
167          if w := opmtch z then return (nil . simp w)
168           else if null cdr y then return (t . mksq(z,1));
169          z := 'sum . (prepsq u . y);
170          let sum_last_attempt_rules!*;
171          w:= opmtch z;
172          rule!-list (list sum_last_attempt_rules!*,nil);
173          return (t . if w then simp w else mksq(z,1))>>;
174%         return (t . if w := opmtch z then simp w else mksq(z,1))>>;
175     z := nil ./ 1;
176  a:
177     if dif < 0 then return (t . z);
178     z := addsq(z,subsq(u,list(v . list('plus,lower,dif))));
179     dif := dif - 1;
180     go to a
181   end;
182
183symbolic procedure sum!-subst(u,x,a);
184   if u = x then a
185    else if atom u then u
186    else sum!-subst(car u, x, a) . sum!-subst(cdr u, x, a);
187
188symbolic procedure sum!-df(u,y);
189   begin scalar w,z,upper,lower,dif;
190      dif := nil;
191      if length(y) = 3 then  <<
192         lower := cadr y;
193         upper := caddr y;
194         dif := addsq(simp!* upper, negsq simp!* lower);
195         if denr dif = 1 then
196            if null numr dif
197               then return simp!* sum!-subst(u, car y, upper)
198               else if fixp numr dif then dif := numr dif
199               else dif := nil
200           else dif := nil;
201         if dif and dif <= 0 then return nil ./ 1 >>;
202     if null dif then <<
203          z := 'sum . (u . y);
204          let sum_last_attempt_rules!*;
205          w:= opmtch z;
206          rule!-list (list sum_last_attempt_rules!*,nil);
207          return if w then simp w else mksq(z,1)>>;
208     z := nil ./ 1;
209  a: if dif < 0 then return z;
210     z := addsq(z,simp!* sum!-subst(u, car y, list('plus,lower,dif)));
211     dif := dif - 1;
212     go to a
213   end;
214
215
216%*********************************************************************
217%       Summation by Gosper's algorithm.
218%********************************************************************;
219
220symbolic procedure sum!-sq(u,v);
221   %Argument U : expression of s-q;
222   %         V : kernel;
223   %value        : expression of sq (result of summation.);
224   begin scalar gn,fn,pn,rn,qn,z,k,x;
225      if null numr u then return nil ./ 1;
226      x := setkorder list v;
227      z := reorder numr u;
228      u := z ./ reorder denr u;
229      if !*trsum then <<
230            prin2t " *** Summation by Gosper's algorithm ***";
231            prin2!* "    A(n) = "; sqprint u;terpri!* t;
232            terpri!* t>>;
233      if domainp z or not (mvar z eq v) or red z then
234            <<pn := 1; z := u>>
235       else <<pn := !*p2f lpow z;
236              z := lc z ./ denr u>>;
237      z := quotsq(z,nsubsq(z,v, - 1));
238      gn := gcdf!*(numr z,denr z);
239      if !*trsum then <<
240         prin2!* "A(n)/A(n-1) = ";sqprint z;terpri!* t;
241         prin2!* "GN = ";xprinf(gn,nil,nil);terpri!* t>>;
242      qn := quotf!*(numr z, gn);
243      rn := quotf!*(denr z, gn);
244      if nonpolyp(qn,v) or nonpolyp(rn,v) then go to fail;
245      if !*trsum then <<
246         prin2!* "Initial qn, rn and pn are "; terpri!* t;
247         prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t;
248         prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t;
249         prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>;
250      k := compress explode '!+j;
251      z := integer!-root(resultant(qn,nsubsf(rn,v,k),v),k);
252      if !*trsum then <<
253         prin2 "Root of resultant(q(n),r(n+j)) are "; prin2t z >>;
254      while z do <<
255            k := car z;
256            gn := gcdf!*(qn,nsubsf(rn,v,k));
257            qn := quotf!*(qn,gn);
258            rn := quotf!*(rn,nsubsf(gn,v, -k));
259            while (k := k - 1)>=0 do
260               pn := multf(pn,nsubsf(gn,v, -k));
261            z := cdr z>>;
262      if !*trsum then <<
263         prin2!* "Shift free  qn, rn and pn are";terpri!* t;
264         prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t;
265         prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t;
266         prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>;
267      qn := nsubsf(qn,v,1);
268      if (k := degree!-bound(pn,addf(qn,rn),addf(qn,negf rn),v)) < 0
269         then go to fail;
270      if !*trsum then <<
271         prin2 "DEGREE BOUND is "; prin2t k >>;
272      if not(fn := solve!-fn(k,pn,qn,rn,v)) then go to fail;
273      if !*trsum then <<
274         prin2!* "FN = ";sqprint fn;terpri!* t >>;
275      u := multsq(multsq(qn ./ 1,fn), multsq(u, 1 ./ pn));
276      z := gcdf!*(numr u, denr u);
277      u := quotf!*(numr u, z) ./ quotf!*(denr u,z);
278      if !*trsum then <<
279            prin2t " *** Gosper's algorithm completed ***";
280            prin2!* "    S(n) = "; sqprint u;terpri!* t;
281            terpri!* t>>;
282      setkorder x;
283      return (reorder numr u ./ reorder denr u);
284  fail:
285      if !*trsum then <<
286            prin2t " *** Gosper's algorithm failed ***";
287            terpri!* t>>;
288      setkorder x;
289      return 'failed
290   end;
291
292%*********************************************************************
293%       integer root isolation
294%********************************************************************;
295
296symbolic procedure integer!-root(u,v);
297% Produce a list of all positive integer root of U;
298   begin scalar x,root,n,w;
299        x := setkorder list v;
300        u := reorder u;
301        if domainp u or not(mvar u eq v) then go to a;
302        u := numr cancel(u ./ lc u);
303        w := u;                                 % get trailing term;
304        while not domainp w and mvar w eq v and cdr w do
305            w := cdr w;
306        if (n := degr(w,v)) > 0 then <<
307            w := lc w;
308            while n > 0 do <<
309               root := 0 . root;
310               n := n - 1>>>>;
311        n := dfactors lowcoef w;                % factor tail coeff.
312        w := (v . 1) . 1;
313        while n do <<
314            if not testdivide(u,v,car n) then <<
315                root := car n . root;
316                u := quotf!*(u, (w . - car n))>>
317             else n := cdr n>>;
318  a:
319        setkorder x;
320        return root
321   end;
322
323
324symbolic procedure lowcoef u;
325   begin scalar lst,m;
326        lst := dcoefl u;
327        m := 0;
328  a:
329        if null lst then return m;
330        m := gcdn(m,car lst);
331        if m = 1 then return 1;
332        lst := cdr lst;
333        go to a
334   end;
335
336
337symbolic procedure dcoefl u;
338   if domainp u then if fixp u then list abs u else nil
339    else nconc(dcoefl lc u , dcoefl red u);
340
341symbolic procedure testdivide(u,v,n);
342% Evaluate U at integer point (V = N);
343   begin scalar x;
344 a:
345       if domainp u or not(mvar u eq v) then return addf(u,x);
346       x := addf(multd(expt(n,ldeg u),lc u),x);
347       if (u := red u) then go to a;
348       return x
349   end;
350
351
352%*********************************************************************
353%********************************************************************;
354
355symbolic procedure degree!-bound(pn,u,v,kern);
356% degree bound for fn;
357%       u: q(n+1) + r(n);
358%       v: q(n+1) - r(n);
359   begin scalar lp,l!+, l!-, x,m,k;
360      x := setkorder list kern;
361      u := reorder u;
362      v := reorder v;
363      pn := reorder pn;
364      l!+ := if u then degr(u,kern) else -1;
365      l!- := if v then degr(v,kern) else -1;
366      lp := if pn then degr(pn,kern) else -1;
367      if l!+ <= l!- then <<k := lp - l!-;go to  a>>;
368      k := lp - l!+ + 1;
369      if l!+ > 0 then u := lc u;
370      if l!- > 0 then v := lc v;
371      if l!+ = l!- + 1 and fixp(m := quotf1(multd(-2,v),u)) then
372        k := max(m,k)
373       else if lp = l!- then k := max(k,0);
374   a:
375      setkorder x;
376      return k
377   end;
378
379%*********************************************************************
380%       calculate polynomial f(n) such that
381%       p(n) - q(n+1)*f(n) + r(n)*f(n-1) = 0;
382%********************************************************************;
383
384symbolic procedure solve!-fn(k,pn,qn,rn,v);
385   begin scalar i,fn,x,y,z,u,w,c,clst,flst;
386      c := makevar('c,0);
387      clst := list c;
388      fn := !*k2f c;
389      i := 0;
390      while (i := i + 1) <= k do <<
391         c := makevar('c,i);
392         clst := c . clst;
393         fn := ((v . i) . !*k2f c) . fn>>;
394      z :=
395       addf(pn,
396            addf(negf multf(qn,fn),
397                 multf(rn,nsubsf(fn,v, - 1))));
398      x := setkorder (v . clst);
399      z := reorder z;
400      c := clst;
401      if !*trsum then <<
402         prin2!* "C Equation is";terpri!* t;
403         xprinf(z,nil,nil);terpri!* t >>;
404  a:
405      if domainp z or
406         domainp (y := if mvar z eq v then lc z else z) then
407           go to fail;
408      w := mvar y;
409      if not(w memq clst) then go to fail;
410      if !*trsum then <<
411         prin2!* "C Equation to solve is ";xprinf(y,nil,nil);terpri!* t;
412         prin2!* " w.r.t ";xprinf(!*k2f w,nil,nil);terpri!* t >>;
413      u := gcdf!*(red y , lc y);
414      u := quotf!*(negf red y, u) ./ quotf!*(lc y, u);
415      flst := (w . u) . flst;
416      z := subst!-cn(z,w,u);
417      if !*trsum then <<
418         xprinf(!*k2f w,nil,nil);prin2!* " := ";sqprint u;terpri!* t >>;
419      clst := deleteq(clst,w);
420      if z then go to a;
421      setkorder c;
422      fn := reorder fn;
423      u := 1;
424      while not domainp fn and mvar fn memq c do <<
425         w := mvar fn;
426         z := atsoc(w,flst);
427         fn := subst!-cn(fn,w,if z then cdr z);
428         if z then u := multf(u,denr cdr z);
429         z := gcdf!*(fn,u);
430         fn := quotf!*(fn,z);
431         u := quotf!*(u,z)>>;
432      setkorder x;
433      return cancel(reorder fn ./ reorder u);
434 fail:
435     if !*trsum then <<
436        prin2t "Fail to solve C equation.";
437        prin2!* "Z := ";xprinf(z,nil,nil);terpri!* t >>;
438     setkorder x;
439     return nil
440   end;
441
442
443symbolic procedure subst!-cn(u,v,x);
444   begin scalar z;
445      z := setkorder list v;
446      u := reorder u;
447      if not domainp u and mvar u eq v then
448         if x then u := addf(multf(lc u,reorder numr x),
449                             multf(red u,reorder denr x))
450          else u := red u;
451      setkorder z;
452      return reorder u
453   end;
454
455
456
457symbolic procedure makevar(id,n);
458  compress nconc(explode id, explode n);
459
460symbolic procedure deleteq(u,x);
461   if null u then nil
462    else if car u eq x then cdr u
463    else car u . deleteq(cdr u, x);
464
465
466symbolic procedure nsubsf(u,kern,i);
467% ARGUMENT U : expression of sf;
468%          KERN : kernel;
469%          I : integer or name of integer variable;
470% value        : expression of sf;
471   begin scalar x,y,z,n;
472      if null i or i = 0 then return u;
473      x := setkorder list kern;
474      u := reorder u;
475      y := addf(!*k2f kern,
476                if fixp i then i else !*k2f i);
477      z := nil;
478   a:
479      if domainp u or not(mvar u eq kern) then goto b;
480      z := addf(z,lc u);
481      n := degr(u,kern) - degr(red u,kern);
482      u := red u;
483   a1:
484      if n <= 0 then goto a;
485      z := multf(z,y);
486      n := n - 1;
487      go to a1;
488   b:
489      z := addf(z,u);
490      setkorder x;
491      return reorder z
492   end;
493
494
495symbolic procedure nsubsq(u,kern,i);
496% ARGUMENT U : expression of sq;
497%       KERN : kernel;
498%          I : integer or name of integer variable;
499% value        : expression of sq;
500   subsq(u,list(kern . list('plus, kern, i)));
501
502
503%*********************************************************************
504%       dependency check
505%********************************************************************;
506
507symbolic procedure nonpolyp(u,v);
508% check U is not a polynomial of V;
509   if domainp u then nil
510    else (not(mvar u eq v) and depend!-p(mvar u,v))
511         or nonpolyp(lc u,v) or nonpolyp(red u,v);
512
513
514symbolic procedure depend!-sq(u,v);
515   depend!-f(numr u,v) or depend!-f(denr u,v);
516
517
518symbolic procedure depend!-f(u,v);
519   if domainp u then nil
520    else depend!-p(mvar u,v) or
521         depend!-f(lc u,v) or
522         depend!-f(red u,v);
523
524
525symbolic procedure depend!-p(u,v);
526   if u eq v then t
527    else if atom u then nil
528    else if not atom car u then depend!-f(u,v)
529    else if car u eq '!*sq then depend!-sq(cadr u, v)
530    else depend!-l(cdr u, v);
531
532symbolic procedure depend!-l(u,v);
533   if null u then nil
534    else if depend!-sq(simp car u, v) then t
535    else depend!-l(cdr u,v);
536
537
538%*********************************************************************
539%       term splitting
540%********************************************************************;
541
542symbolic procedure sum!-term!-split(u,v);
543   begin scalar y,z,klst,lst,x;
544      x := setkorder list v;
545      z := qremf(reorder numr u, y := reorder denr u);
546      klst := kern!-list(car z,v);
547      lst := termlst(car z, 1 ./ 1, klst);
548      klst := kern!-list(cdr z,v);
549      if depend!-f(y,v) then klst := deleteq(klst,v);
550      lst := append(lst, termlst(cdr z, 1 ./ y, klst));
551      setkorder x;
552      return lst
553   end;
554
555
556symbolic procedure kern!-list(u,v);
557% Returns list of kernels that depend on V;
558   begin scalar x;
559      for each j in kernels u do if depend!-p(j,v) then x := j . x;
560      return x
561   end;
562
563symbolic procedure termlst(u,v,klst);
564   begin scalar x,kern,lst;
565      if null u then return nil
566       else if null klst or domainp u
567        % Preserve order for noncom.
568        then return list multsq(v,!*f2q u);
569      kern := car klst;
570      klst := cdr klst;
571      x := setkorder list kern;
572      u := reorder u;
573      v := reorder(numr v) ./ reorder(denr v);
574      while not domainp u and mvar u eq kern do <<
575         lst := nconc(termlst(lc u, multsq(!*p2q lpow u, v),klst),lst);
576         u := red u>>;
577      if u then lst := nconc(termlst(u,v,klst),lst);
578      setkorder x;
579      return lst
580   end;
581
582
583%*********************************************************************
584%       Express trigonometric functions (such as sin, cos ..)
585%       by exponentials.
586%********************************************************************;
587
588symbolic procedure trig!-to!-expsq(u,v);
589   multsq(trig!-to!-expf(numr u,v),
590          invsq trig!-to!-expf(denr u,v));
591
592
593symbolic procedure trig!-to!-expf(u,v);
594   if domainp u then u ./ 1
595    else addsq(multsq(trig!-to!-expp(lpow u,v),
596                      trig!-to!-expf(lc u,v)),
597               trig!-to!-expf(red u,v));
598
599
600symbolic procedure trig!-to!-expp(u,v);
601   begin scalar !*combineexpt,w,x,z,n,wi;
602      % We don't want to combine expt terms here, since the code
603      % depends on the terms being separate.
604      n := cdr u;          % integer power;
605      z := car u;          % main variable;
606      if atom z or not atom (x := car z) or not depend!-p(z,v) then
607           return !*p2q u;
608      if x memq '(sin cos tan sec cosec cot) then <<
609          !*trig!-to!-exp := t;
610          w := multsq(!*k2q 'i, simp!* cadr z);
611          w := simp!* list('expt,'e, mk!*sq w);
612%         W := SIMP LIST('EXPT,'E, 'TIMES . ( 'I . CDR Z));
613          wi := invsq w;
614          if x eq 'sin then
615              w := multsq(addsq(w ,negsq wi),
616                          1 ./ list(('i .** 1) .* 2))
617           else if x eq 'cos then
618              w := multsq(addsq(w, wi), 1 ./ 2)
619           else if x eq 'tan then
620              w := multsq(addsq(w,negsq wi),
621                          invsq addsq(w,wi))
622           else if x eq 'sec then
623              w := multsq(2 ./ 1, invsq addsq(w, wi))
624           else if x eq 'cosec then
625              w := multsq(list(('i .** 1) .* 2),
626                          invsq addsq(w, negsq wi))
627           else
628              w := multsq(addsq(w, wi),
629                          invsq addsq(w, negsq wi))
630                     >>
631       else if x memq '(sinh cosh tanh sech cosech coth) then <<
632          !*trig!-to!-exp := t;
633          w := simp!* list('expt,'e,cadr z);
634          wi := invsq w;
635          if x eq 'sinh then
636              w := multsq(addsq(w,negsq wi), 1 ./ 2)
637           else if x eq 'cosh then
638              w := multsq(addsq(w,wi), 1 ./ 2)
639           else if x eq 'tanh then
640              w := multsq(addsq(w,negsq wi),
641                          invsq addsq(w,wi))
642           else if x eq 'sech then
643              w := multsq(2 ./ 1, invsq addsq(w, wi))
644           else if x eq 'cosech then
645              w := multsq(2 ./ 1, invsq addsq(w, negsq wi))
646           else
647              w := multsq(addsq(w,wi),
648                          invsq addsq(w, negsq wi))
649                 >>
650       else return !*p2q u;
651      return exptsq(w,n)
652   end;
653
654%*********************************************************************
655%       Inverse of trig!-to!-exp.
656%       Express exponentials in terms of trigonometric functions
657%                (sin, cos, sinh and cosh)
658%                               Wed Dec. 17, 1986 by F. Kako;
659%********************************************************************;
660
661symbolic procedure exp!-to!-trigsq u;
662   multsq(exp!-to!-trigf numr u,
663          invsq exp!-to!-trigf denr u);
664
665symbolic procedure exp!-to!-trigf u;
666   begin scalar v,v1,x,y,n;
667      u := termlst1(u,1,nil ./1);
668      v := nil;
669   a:
670      if null u then go to b;
671      x := caar u;
672      y := cdar u;
673      u := cdr u;
674   a1:
675      if u and y = cdar u then <<
676         x := addf(x,caar u);
677         u := cdr u;
678         go to a1>>;
679      v := (x . y) . v;
680      go to a;
681   b:
682      v1 := reverse v;
683      n := length v;
684      u := nil ./ 1;
685   c:
686      if n = 0 then return u
687       else if n = 1 then
688          return addsq(u,
689                       multsq(!*f2q caar v,
690                              simp!* list('expt,'e,mk!*sq cdar v)));
691      u := addsq(u,exp!-to!-trigl(caar v1,caar v,cdar v1,cdar v));
692      v := cdr v;
693      v1 := cdr v1;
694      n := n - 2;
695      go to c
696   end;
697
698symbolic procedure exp!-to!-trigl(a,b,c,d);
699%       A*E**C + B*E**D
700%               -->
701%       ((A+B)*COSH((C-D)/2)+(A-B)*SINH((C-D)/2))*E**((C+D)/2);
702%       A, B: sf;
703%       C, D: sq;
704   begin scalar x,y,z;
705      x := !*f2q addf(a,b);
706      y := !*f2q addf(a, negf b);
707      z := multsq(addsq(c,negsq d), 1 ./ 2);
708      z := real!-imag!-sincos z;
709      return multsq(simp!* list('expt,'e,
710                                mk!*sq multsq(addsq(c,d),1 ./ 2)),
711                    addsq(multsq(x, cdr z),
712                          multsq(y, car z)))
713   end;
714
715
716symbolic procedure termlst1(u,v,w);
717%ARGUMENT U : sf;
718%         V : sf;
719%         W : sq;
720%value      : list of (sf . sq);
721   begin scalar x,y;
722      if null u then return nil
723       else if domainp u then return list (multf(u,v) . w);
724      x := mvar u;
725      y := if atom x or not(car x eq 'expt) or not(cadr x eq 'e) then
726              termlst1(lc u,multf(!*p2f lpow u,v),w)
727            else termlst1(lc u,v,addsq(w,multsq(simp!* caddr x,
728                                                ldeg u ./ 1)));
729      return nconc(y,termlst1(red u,v,w))
730   end;
731
732
733% These can be found in Abramowitz-Stegun (27.8.6 Summable Series), and
734% were suggested by Winfried Neun.
735
736algebraic;
737
738let {sum(sin(~n*~tt)/n,~n,1,infinity) => (pi - tt)/2,
739     sum(sin(~n*~tt)/(~n)^3,~n,1,infinity) =>
740             pi^2*tt/6 - pi*tt^2/4 + tt^3/12,
741     sum(sin(~n*~tt)/(~n)^5,~n,1,infinity) =>
742             pi^4*tt/90 - pi^2*tt^3/36 + pi*tt^4/48-tt^5/240}$
743
744let {sum(cos(~n*~tt)/(~n),~n,1,infinity) =>  -log(2*sin(tt/2)),
745     sum(cos(~n*~tt)/(~n)^2,~n,1,infinity) =>
746             pi^2/6 - pi*tt/2 + tt^2/4,
747     sum(cos(~n*~tt)/(~n)^4,~n,1,infinity) =>
748             pi^4/90 - pi^2*tt^2/12 + pi*tt^3/12-tt^4/48}$
749
750endmodule;
751
752end;
753