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