1module tayutils; 2 3% Redistribution and use in source and binary forms, with or without 4% modification, are permitted provided that the following conditions are met: 5% 6% * Redistributions of source code must retain the relevant copyright 7% notice, this list of conditions and the following disclaimer. 8% * Redistributions in binary form must reproduce the above copyright 9% notice, this list of conditions and the following disclaimer in the 10% documentation and/or other materials provided with the distribution. 11% 12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 16% CONTRIBUTORS 17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 23% POSSIBILITY OF SUCH DAMAGE. 24% 25 26 27%***************************************************************** 28% 29% Utility functions that operate on Taylor kernels 30% 31%***************************************************************** 32 33 34exports add!-degrees, add!.comp!.tp!., check!-for!-cst!-taylor, 35 comp!.tp!.!-p, delete!-superfluous!-coeffs, enter!-sorted, 36 exceeds!-order, exceeds!-order!-variant, find!-non!-zero, 37 get!-cst!-coeff, inv!.tp!., is!-neg!-pl, min2!-order, 38 mult!.comp!.tp!., rat!-kern!-pow, replace!-next, 39 subtr!-degrees, subtr!-tp!-order, taydegree!<, 40 taydegree!-strict!<!=, taymincoeff, tayminpowerlist, 41 taylor!*!-constantp, taylor!*!-nzconstantp, taylor!*!-onep, 42 taylor!*!-zerop, taytpmin2, tp!-greaterp, truncate!-coefflist, 43 truncate!-taylor!*; 44 45imports 46 47% from the REDUCE kernel: 48 ./, gcdn, geq, lastpair, mkrn, neq, nth, numr, reversip, 49 50% from the header module: 51 !*tay2q, get!-degree, get!-degreelist, make!-cst!-powerlist, 52 make!-taylor!*, nzerolist, prune!-coefflist, taycfpl, taycfsq, 53 taycoefflist, taygetcoeff, tayflags, taylor!:, tayorig, 54 taytemplate, taytpelnext, taytpelorder, taytpelpoint, 55 taytpelvars, tpdegreelist, tpnextlist, 56 57% from module Tayintro: 58 confusion, 59 60% from module TayUtils: 61 taycoefflist!-zerop; 62 63 64 65symbolic procedure add!-degrees (dl1, dl2); 66 % 67 % calculates the element-wise sum of two degree lists dl1, dl2. 68 % 69 taylor!: 70 begin scalar dl,u,v,w; 71 while dl1 do << 72 u := car dl1; 73 v := car dl2; 74 w := nil; 75 while u do << 76 w := (car u + car v) . w; 77 u := cdr u; 78 v := cdr v>>; 79 dl := reversip w . dl; 80 dl1 := cdr dl1; 81 dl2 := cdr dl2>>; 82 return reversip dl 83 end; 84 85symbolic procedure subtr!-degrees(dl1,dl2); 86 % 87 % calculates the element-wise difference of two degree lists dl1, dl2. 88 % 89 taylor!: 90 begin scalar dl,u,v,w; 91 while dl1 do << 92 u := car dl1; 93 v := car dl2; 94 w := nil; 95 while u do << 96 w := (car u - car v) . w; 97 u := cdr u; 98 v := cdr v>>; 99 dl := reversip w . dl; 100 dl1 := cdr dl1; 101 dl2 := cdr dl2>>; 102 return reversip dl 103 end; 104 105symbolic procedure find!-non!-zero pl; 106 % 107 % pl is a power list. Returns pair (n . m) of position of first 108 % non zero degree. 109 % 110 begin scalar u; integer n, m; 111 n := 1; 112 loop: 113 m := 1; 114 u := car pl; 115 loop2: 116 if not (car u = 0) then return (n . m); 117 u := cdr u; 118 m := m + 1; 119 if not null u then goto loop2; 120 pl := cdr pl; 121 n := n + 1; 122 if null pl then confusion 'find!-non!-zero; 123 goto loop 124 end$ 125 126 127symbolic procedure lcmn(n,m); 128 % 129 % returns the least common multiple of two integers m,n. 130 % 131 n*(m/gcdn(n,m)); 132 133symbolic inline procedure get!-denom expo; 134 if atom expo then 1 else cddr expo; 135 136symbolic procedure get!-denom!-l expol; 137 <<for each el in cdr expol do 138 result := lcmn(result,get!-denom el); 139 result>> 140 where result := get!-denom car expol; 141 142symbolic procedure get!-denom!-ll(dl,pl); 143 if null dl then nil 144 else lcmn(car dl,get!-denom!-l car pl) 145 . get!-denom!-ll(cdr dl, cdr pl); 146 147symbolic procedure smallest!-increment cfl; 148 if null cfl then confusion 'smallest!-increment 149 else begin scalar result; 150 result := for each l in taycfpl car cfl collect get!-denom!-l l; 151 for each el in cdr cfl do 152 result := get!-denom!-ll(result,taycfpl el); 153 return for each n in result collect if n=1 then n else mkrn(1,n); 154 end; 155 156 157symbolic procedure common!-increment(dl1,dl2); 158 begin scalar result,l; 159 loop: 160 l := lcmn(get!-denom car dl1,get!-denom car dl2); 161 result := (if l=1 then l else mkrn(1,l)) . result; 162 dl1 := cdr dl1; 163 dl2 := cdr dl2; 164 if not null dl1 then goto loop 165 else if not null dl2 then confusion 'common!-increment 166 else return reversip result 167 end; 168 169symbolic procedure min2!-order(nextlis,ordlis,dl); 170 % 171 % (List of Integers, List of Integers, TayPowerList) -> Boolean 172 % 173 % nextlis is the list of TayTpElNext numbers, 174 % ordlis the list of TayTpElOrder numbers, 175 % dl the degreelist of a coefficient. 176 % Dcecrease the TayTpElNext number if the degree is greater than 177 % the order, but smaller than the next. 178 % Returns the corrected nextlis. 179 % 180 if null nextlis then nil 181 else (taylor!: (if dg > car ordlis then min2(dg,car nextlis) 182 else car nextlis) where dg := get!-degree car dl) 183 . min2!-order(cdr nextlis,cdr ordlis,cdr dl); 184 185symbolic procedure replace!-next(tp,nl); 186 % 187 % Given a template and a list of exponents, returns a template 188 % with the next part replaced. 189 % 190 if null tp then nil 191 else {taytpelvars car tp, 192 taytpelpoint car tp, 193 taytpelorder car tp, 194 car nl} 195 . replace!-next(cdr tp,cdr nl); 196 197symbolic procedure comp!.tp!.!-p (u, v); 198 % 199 % Checks templates of Taylor kernels u and v for compatibility, 200 % i.e. whether variables and expansion points match. 201 % Returns t if possible. 202 % 203 begin; 204 u := taytemplate u; v := taytemplate v; 205 if length u neq length v then return nil; 206 loop: 207 if not (taytpelvars car u = taytpelvars car v 208 and taytpelpoint car u = taytpelpoint car v) 209 then return nil; 210 u := cdr u; v := cdr v; 211 if null u then return t; 212 goto loop 213 end$ 214 215symbolic procedure add!.comp!.tp!.(u,v); 216 % 217 % Checks templates of Taylor kernels u and v for compatibility 218 % when adding them, i.e. whether variables and expansion points 219 % match. 220 % Returns either a list containing a new Taylor template whose 221 % degrees are the minimum of the corresponding degrees of u and v, 222 % or nil if variables or expansion point(s) do not match. 223 % 224 taylor!: 225 begin scalar w; 226 u := taytemplate u; v := taytemplate v; 227 if length u neq length v then return nil; 228 if null u then return {nil}; 229 loop: 230 if not (taytpelvars car u = taytpelvars car v 231 and taytpelpoint car u = taytpelpoint car v) 232 then return nil 233 else w := {taytpelvars car u, 234 taytpelpoint car u, 235 min2(taytpelorder car u,taytpelorder car v), 236 min2(taytpelnext car u,taytpelnext car v)} 237 . w; 238 u := cdr u; v := cdr v; 239 if null u then return {reversip w}; 240 goto loop 241 end; 242 243symbolic procedure taymindegreel(pl,dl); 244 taylor!: 245 if null pl then nil 246 else min2(get!-degree car pl,car dl) 247 . taymindegreel(cdr pl,cdr dl); 248 249symbolic procedure get!-min!-degreelist cfl; 250 taylor!: 251 if null cfl then confusion 'get!-min!-degreelist 252 else if null cdr cfl then get!-degreelist taycfpl car cfl 253 else taymindegreel(taycfpl car cfl, 254 get!-min!-degreelist cdr cfl); 255 256symbolic procedure mult!.comp!.tp!.(u,v,div!?); 257 % 258 % Checks templates of Taylor kernels u and v for compatibility 259 % when multiplying or dividing them, i.e., whether variables and 260 % expansion points match. The difference to addition is that 261 % in addition to the new template it returns two degreelists 262 % and two nextlists to be used by truncate!-coefflist which 263 % are made up so that the kernels have the same number of terms. 264 % 265 taylor!: 266 begin scalar cf1,cf2,next1,next2,ord1,ord2,w, 267 !#terms!-1,!#terms!-next,dl1,dl2,mindg; 268 cf1 := prune!-coefflist taycoefflist u; 269 if null cf1 then dl1 := nzerolist length taytemplate u 270 else dl1 := get!-min!-degreelist cf1; 271 cf2 := prune!-coefflist taycoefflist v; 272 if null cf2 then dl2 := nzerolist length taytemplate v 273 else dl2 := get!-min!-degreelist cf2; 274 u := taytemplate u; v := taytemplate v; 275 if length u neq length v then return nil; 276 if null u then return {nil,nil,nil,nil,nil}; 277 loop: 278 if not (taytpelvars car u = taytpelvars car v 279 and taytpelpoint car u = taytpelpoint car v) 280 then return nil; 281 mindg := if div!? then car dl1 - car dl2 else car dl1 + car dl2; 282 !#terms!-1 := min2(taytpelorder car u - car dl1, 283 taytpelorder car v - car dl2); 284 !#terms!-next := min2(taytpelnext car u - car dl1, 285 taytpelnext car v - car dl2); 286 ord1 := (!#terms!-1 + car dl1) . ord1; 287 ord2 := (!#terms!-1 + car dl2) . ord2; 288 next1 := (!#terms!-next + car dl1) . next1; 289 next2 := (!#terms!-next + car dl2) . next2; 290 w := {taytpelvars car u,taytpelpoint car u, 291 mindg + !#terms!-1,mindg + !#terms!-next} 292 . w; 293 u := cdr u; v := cdr v; dl1 := cdr dl1; dl2 := cdr dl2; 294 if null u then return {reversip w, 295 reversip ord1, 296 reversip ord2, 297 reversip next1, 298 reversip next2}; 299 goto loop 300 end; 301 302symbolic procedure inv!.tp!. u; 303 % 304 % Checks template of Taylor kernel u for inversion. It returns a 305 % template (to be used by truncate!-coefflist) 306 % which is made up so that the resulting kernel has the correct 307 % number of terms. 308 % 309 taylor!: 310 begin scalar w,cf,!#terms!-1,!#terms!-next,dl,mindg; 311 cf := prune!-coefflist taycoefflist u; 312 if null cf then dl := nzerolist length taytemplate u 313 else dl := get!-degreelist taycfpl car cf; 314 u := taytemplate u; 315 if null u then return {nil,nil}; 316 loop: 317 mindg := - car dl; 318 !#terms!-1 := taytpelorder car u - car dl; 319 !#terms!-next := taytpelnext car u - car dl; 320 w := {taytpelvars car u,taytpelpoint car u,mindg + !#terms!-1, 321 mindg + !#terms!-next} 322 . w; 323 u := cdr u; dl := cdr dl; 324 if null u then return {reversip w}; 325 goto loop 326 end; 327 328symbolic inline procedure taycoeff!-before(cc1,cc2); 329 % 330 % (TayCoeff, TayCoeff) -> Boolean 331 % 332 % returns t if coeff cc1 is ordered before cc2 333 % both are of the form (degreelist . sq) 334 % 335 taydegree!<(taycfpl cc1,taycfpl cc2); 336 337symbolic procedure taydegree!<(u,v); 338 % 339 % (TayPowerList, TayPowerList) -> Boolean 340 % 341 % returns t if coefflist u is ordered before v 342 % 343 taylor!: 344 begin scalar u1,v1; 345 loop: 346 u1 := car u; 347 v1 := car v; 348 loop2: 349 if car u1 > car v1 then return nil 350 else if car u1 < car v1 then return t; 351 u1 := cdr u1; 352 v1 := cdr v1; 353 if not null u1 then go to loop2; 354 u := cdr u; 355 v := cdr v; 356 if not null u then go to loop 357 end; 358 359symbolic procedure taydegree!-strict!<!=(u,v); 360 % 361 % (TayPowerList, TayPowerList) -> Boolean 362 % 363 % returns t if every component coefflist u is less or equal than 364 % respective component of v 365 % 366 taylor!: 367 begin scalar u1,v1; 368 loop: 369 u1 := car u; 370 v1 := car v; 371 loop2: 372 if car u1 > car v1 then return nil; 373 u1 := cdr u1; 374 v1 := cdr v1; 375 if not null u1 then go to loop2; 376 u := cdr u; 377 v := cdr v; 378 if not null u then go to loop; 379 return t 380 end; 381 382symbolic procedure exceeds!-order(ordlis,cf); 383 % 384 % (List of Integers, TayPowerlist) -> Boolean 385 % 386 % Returns t if the degrees in coefficient cf are greater or 387 % equal than those in the degreelist ordlis 388 % 389 if null ordlis then nil 390 else taylor!:(get!-degree car cf >= car ordlis) 391 or exceeds!-order(cdr ordlis,cdr cf); 392 393symbolic procedure exceeds!-order!-variant(ordlis,cf); 394 % 395 % (List of Integers, TayPowerlist) -> Boolean 396 % 397 % Returns t if the degrees in coefficient cf are greater or 398 % equal than those in the degreelist ordlis 399 % 400 if null ordlis then nil 401 else taylor!:(get!-degree car cf > car ordlis) 402 or exceeds!-order!-variant(cdr ordlis,cdr cf); 403 404symbolic procedure enter!-sorted (u, alist); 405 % 406 % (TayCoeff, TayCoeffList) -> TayCoeffList 407 % 408 % enters u into the alist alist according to the standard 409 % ordering for the car part 410 % 411 if null alist then {u} 412 else if taycoeff!-before (u, car alist) then u . alist 413 else car alist . enter!-sorted (u, cdr alist)$ 414 415symbolic procedure delete!-superfluous!-coeffs(cflis,pos,n); 416 % 417 % (TayCoeffList, Integer, Integer) -> TayCoeffList 418 % 419 % This procedure deletes all coefficients of a TayCoeffList cflis 420 % whose degree in position pos exceeds n. 421 % 422 taylor!: 423 for each cc in cflis join 424 (if get!-degree nth(taycfpl cc,pos) > n then nil else {cc}); 425 426symbolic procedure truncate!-coefflist (cflis, dl); 427 % 428 % (TayCoeffList, List of Integers) -> TayCoeffList 429 % 430 % Deletes all coeffs from coefflist cflis that are equal or greater 431 % in degree than the corresponding degree in the degreelist dl. 432 % 433 begin scalar l; 434 for each cf in cflis do 435 if not exceeds!-order (dl, taycfpl cf) then l := cf . l; 436 return reversip l 437 end; 438 439symbolic procedure taytp!-min2(tp1,tp2); 440 % 441 % finds minimum (w.r.t. Order and Next parts) of compatible 442 % templates tp1 and tp2 443 % 444 taylor!: 445 if null tp1 then nil 446 else if not (taytpelvars car tp1 = taytpelvars car tp2 and 447 taytpelpoint car tp1 = taytpelpoint car tp2) 448 then confusion 'taytpmin2 449 else {taytpelvars car tp1,taytpelpoint car tp2, 450 min2(taytpelorder car tp1,taytpelorder car tp2), 451 min2(taytpelnext car tp1,taytpelnext car tp2)} 452 . taytp!-min2(cdr tp1,cdr tp2); 453 454 455symbolic procedure truncate!-taylor!*(tay,ntp); 456 % 457 % tcl is a coefflist for template otp 458 % truncate it to coefflist for template ntp 459 % 460 taylor!: 461 begin scalar nl,ol,l,tp,tcl,otp; 462 tcl := taycoefflist tay; 463 otp := taytemplate tay; 464 tp := for each pp in pair(ntp,otp) collect 465 {taytpelvars car pp, 466 taytpelpoint car pp, 467 min2(taytpelorder car pp,taytpelorder cdr pp), 468 min2(taytpelnext car pp,taytpelnext cdr pp)}; 469 nl := tpnextlist tp; 470 ol := tpdegreelist tp; 471 for each cf in tcl do 472 if not null numr taycfsq cf 473% then ((if not exceeds!-order(nl,pl) then l := cf . l 474% else nl := min2!-order(nl,ol,pl)) 475 then ((if not exceeds!-order!-variant(ol,pl) then l := cf . l 476 else nl := min2!-order(nl,ol,pl)) 477 where pl := taycfpl cf); 478 return make!-taylor!*(reversip l,replace!-next(tp,nl), 479 tayorig tay,tayflags tay) 480 end; 481 482symbolic procedure tp!-greaterp(tp1,tp2); 483 % 484 % Given two templates tp1 and tp2 with matching variables and 485 % expansion points this function returns t if the expansion 486 % order wrt at least one variable is greater in tp1 than in tp2. 487 % 488 if null tp1 then nil 489 else taylor!: (taytpelorder car tp1 > taytpelorder car tp2) 490 or tp!-greaterp(cdr tp1,cdr tp2); 491 492symbolic procedure subtr!-tp!-order(tp1,tp2); 493 % 494 % Given two templates tp1 and tp2 with matching variables and 495 % expansion points this function returns the difference in their 496 % orders. 497 % 498 taylor!: 499 if null tp1 then nil 500 else (taytpelorder car tp1 - taytpelorder car tp2) 501 . subtr!-tp!-order(cdr tp1,cdr tp2); 502 503 504COMMENT Procedures to non-destructively modify Taylor templates; 505 506symbolic procedure addto!-all!-taytpelorders(tp,nl); 507 taylor!: 508 if null tp then nil 509 else {taytpelvars car tp, 510 taytpelpoint car tp, 511 taytpelorder car tp + car nl, 512 taytpelnext car tp + car nl} . 513 addto!-all!-taytpelorders(cdr tp,cdr nl); 514 515symbolic procedure taymincoeff cflis; 516 % 517 % Returns degree of first non-zero coefficient 518 % or 0 if there isn't any. 519 % 520 if null cflis then 0 521 else if not null numr taycfsq car cflis 522 then get!-degree car taycfpl car cflis 523 else taymincoeff cdr cflis; 524 525symbolic procedure tayminpowerlist cflis; 526 % 527 % Returns degreelist of first non-zero coefficient of TayCoeffList 528 % cflis or a list of zeroes if there isn't any. 529 % 530 if null cflis then confusion 'tayminpowerlist 531 else tayminpowerlist1(cflis,length taycfpl car cflis); 532 533symbolic procedure tayminpowerlist1(cflis,l); 534 if null cflis then nzerolist l 535 else if null numr taycfsq car cflis 536 then tayminpowerlist1(cdr cflis,l) 537 else get!-degreelist taycfpl car cflis; 538 539symbolic procedure get!-cst!-coeff tay; 540 taygetcoeff(make!-cst!-powerlist taytemplate tay,taycoefflist tay); 541 542symbolic procedure taylor!*!-constantp tay; 543 taylor!*!-constantp1(make!-cst!-powerlist taytemplate tay, 544 taycoefflist tay); 545 546symbolic procedure taylor!*!-constantp1(pl,tcf); 547 if null tcf then t 548 else if taycfpl car tcf = pl 549 then taycoefflist!-zerop cdr tcf 550 else if not null numr taycfsq car tcf then nil 551 else taylor!*!-constantp1(pl,cdr tcf); 552 553symbolic procedure check!-for!-cst!-taylor tay; 554 begin scalar pl,tc; 555 pl := make!-cst!-powerlist taytemplate tay; 556 tc := taycoefflist tay; 557 return if taylor!*!-constantp1(pl,tc) then taygetcoeff(pl,tc) 558 else !*tay2q tay 559 end; 560 561symbolic procedure taylor!*!-nzconstantp tay; 562 taylor!*!-nzconstantp1(make!-cst!-powerlist taytemplate tay, 563 taycoefflist tay); 564 565symbolic procedure taylor!*!-nzconstantp1(pl,tcf); 566 if null tcf then nil 567 else if taycfpl car tcf = pl 568 then if null numr taycfsq car tcf then nil 569 else taycoefflist!-zerop cdr tcf 570 else if taycfpl car tcf neq pl and 571 not null numr taycfsq car tcf 572 then nil 573 else taylor!*!-nzconstantp1(pl,cdr tcf); 574 575symbolic procedure taylor!*!-onep tay; 576 taylor!-onep1(make!-cst!-powerlist taytemplate tay,taycoefflist tay); 577 578symbolic procedure taylor!-onep1(pl,tcf); 579 if null tcf then nil 580 else if taycfpl car tcf = pl 581 then if taycfsq car tcf = (1 ./ 1) 582 then taycoefflist!-zerop cdr tcf 583 else nil 584 else if null numr taycfsq car tcf 585 then taylor!*!-nzconstantp1(pl,cdr tcf) 586 else nil; 587 588symbolic procedure is!-neg!-pl pl; 589 % 590 % Returns t if any of the exponents in pl is negative. 591 % 592 taylor!: 593 if null pl then nil 594 else if get!-degree car pl < 0 then t 595 else is!-neg!-pl cdr pl; 596 597symbolic procedure rat!-kern!-pow(x,pos); 598 % 599 % check that s.f. x is a kernel to a rational power. 600 % if pos is t allow only positive exponents. 601 % returns pair (kernel . power) 602 % 603 begin scalar y; integer n; 604 if domainp x or not null red x or not (lc x=1) then return nil; 605 n := ldeg x; 606 x := mvar x; 607 taylor!: 608 if eqcar(x,'sqrt) then return (cadr x . mkrn(1,2)*n) 609 else if eqcar(x,'expt) and (y := simprn{caddr x}) 610 then if null pos or (y := car y)>0 611 then return (cadr x . (y*n)) 612 else return nil 613 else return (x . n) 614 end; 615 616endmodule; 617 618end; 619