1module tayintrf; 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% The interface to the REDUCE simplificator 30% 31%***************************************************************** 32 33 34exports simptaylor, simptaylor!*, taylorexpand$ 35 36imports 37 38% from the REDUCE kernel: 39 !*f2q, aconc!*, denr, depends, diffsq, eqcar, kernp, lastpair, 40 leq, lprim, mkquote, mksq, multsq, mvar, neq, nth, numr, over, 41 prepsq, revlis, reversip, simp, simp!*, subs2, subsq, typerr, 42 43% from the header module: 44 !*tay2q, get!-degree, has!-taylor!*, has!-tayvars, 45 make!-taylor!*, multintocoefflist, resimptaylor, taycfpl, 46 taycfsq, taycoefflist, tayflags, taymakecoeff, tayorig, 47 taytemplate, taytpelorder, taytpelpoint, 48 taylor!-kernel!-sq!-p, taymincoeff, 49 50% from module Tayintro: 51 replace!-nth, taylor!-error, var!-is!-nth, 52 53% from module TayExpnd: 54 taylorexpand, 55 56% from module Tayutils: 57 delete!-superfluous!-coeffs, 58 59% from module taybasic: 60 invtaylor1, quottaylor1, 61 62% from module Tayconv: 63 preptaylor!*; 64 65 66fluid '(!*backtrace !*precise !*tayinternal!* !*taylorkeeporiginal !*taylorautocombine taynomul!* 67 frlis!* subfg!*); 68 69global '(mul!*); 70 71COMMENT The following statement forces all expressions to be 72 re-simplified if the switch `taylorautocombine' is set to on, 73 unfortunately, this is not always sufficient; 74 75put ('taylorautocombine, 'simpfg, '((t (rmsubs)))); 76 77 78COMMENT Interface to the fkern mechanism that makes kernels unique and 79 stores them in the klist property; 80 81symbolic procedure tayfkern u; 82 begin scalar x,y; 83 if !*tayinternal!* then return u; 84 % rest of code borrowed from fkern 85 y := get('taylor!*,'klist); 86% Here I will do the potentially slow search using assoc rather than 87% hash table lookup. 88 x := assoc(u,y); 89 if null x 90 then <<x := list(u,nil); 91#if (memq 'csl lispsystem!*) 92 puthash(u, kernhash, x); 93#endif 94 y := ordad(x,y); 95 kprops!* := union('(taylor!*),kprops!*); 96 put('taylor!*,'klist,y)>>; 97 return x 98 end; 99 100put('taylor!*, 'fkernfn, 'tayfkern); 101 102symbolic procedure taysimpsq!-from!-mul u; 103 (taysimpsq u) where taynomul!* := t; 104 105symbolic procedure simptaylor u; 106 % 107 % (PrefixForm) -> s.q. 108 % 109 % This procedure is called directly by the simplifier. 110 % Its argument list must be of the form 111 % (exp, [var, var0, deg, ...]), 112 % where [...] indicate one or more occurences. 113 % This means that exp is to be expanded w.r.t var about var0 114 % up to degree deg. 115 % var may also be a list of variables, which means that the 116 % the expansion takes place in a homogeneous way. 117 % If var0 is the special atom infinity var is replaced by 1/var 118 % and the result expanded about 0. 119 % 120 % This procedure returns the input if it cannot expand the expression. 121 % 122 if remainder(length u,3) neq 1 123 then taylor!-error('wrong!-no!-args,'taylor) 124 else if null subfg!* then mksq('taylor . u,1) 125 else begin scalar !*precise,arglist,degree,f,ll,result,var,var0; 126 % 127 % Allow automatic combination of Taylor kernels. 128 % 129 if !*taylorautocombine and not taynomul!* and not ('taysimpsq memq mul!*) 130 then mul!* := aconc!*(mul!*,'taysimpsq!-from!-mul); 131 % 132 % First of all, call the simplifier on exp (i.e. car u), 133 % 134 f := simp!* car u; 135 u := revlis cdr u; % reval instead of simp!* to handle lists 136 arglist := u; 137 % 138 % then scan the rest of the argument list. 139 % 140 while not null arglist do 141 << var := car arglist; 142 var := if eqcar(var,'list) then cdr var else {var}; 143 % In principle one should use !*a2k 144 % but typerr (maprin) does not correctly print the atom nil 145 for each el in var collect begin 146 el := simp!* el; 147 if kernp el then return mvar numr el 148 else typerr(prepsq el,'kernel) 149 end; 150 var0 := cadr arglist; 151 degree := caddr arglist; 152 if not fixp degree 153 then typerr(degree,"order of Taylor expansion"); 154 arglist := cdddr arglist; 155 ll := {var,var0,degree,degree + 1} . ll>>; 156 % 157 % Now ll is a Taylor template, i.e. of the form 158 % ((var_1 var0_1 deg1 next_1) (var_2 var0_2 deg_2 next_2) ...) 159 % 160 result := taylorexpand(f,reversip ll); 161 % 162 % If the result does not contain a Taylor kernel, return the input. 163 % 164 return if has!-taylor!* result then result 165 else mksq('taylor . prepsq f . u,1) 166 end; 167 168put('taylor,'simpfn,'simptaylor)$ 169 170 171%symbolic procedure taylorexpand (f, ll); 172% % 173% % f is a s.q. that is expanded according to the list ll 174% % which has the form 175% % ((var_1 var0_1 deg1) (var_2 var0_2 deg_2) ...) 176% % 177% begin scalar result; 178% result := f; 179% for each el in ll do 180% % 181% % taylor1 is the function that does the real work 182% % 183% result := !*tay2q taylor1 (result, car el, cadr el, caddr el); 184% if !*taylorkeeporiginal then set!-TayOrig (mvar numr result, f); 185% return result 186% end$ 187 188 189symbolic procedure taylor1 (f, varlis, var0, n); 190 % 191 % Taylor expands s.q. f w.r.t. varlis about var0 up to degree n. 192 % var is a list of kernels, which means that the expansion 193 % takes place in a homogeneous way if there is more than one 194 % kernel. 195 % If var0 is the special atom infinity all kernels in varlis are 196 % replaced by 1/kernel. The result is then expanded about 0. 197 % 198 taylor1sq (if var0 eq 'infinity 199 then subsq (f, 200 for each krnl in varlis collect 201 (krnl . list ('quotient, 1, krnl))) 202 else f, 203 varlis, var0, n)$ 204 205symbolic procedure taylor1sq (f, varlis, var0, n); 206 % 207 % f is a standard quotient, value is a Taylor kernel. 208 % 209 % First check if it is a Taylor kernel 210 % 211 if taylor!-kernel!-sq!-p f 212 then if has!-tayvars(mvar numr f,varlis) 213 % 214 % special case: f has already been expanded w.r.t. var. 215 % 216 then taylorsamevar (mvar numr f, varlis, var0, n) 217 else begin scalar y, z; 218 f := mvar numr f; 219 % 220 % taylor2 returns a list of the form 221 % ((deg1 . coeff1) (deg2 . coeff2) ... ) 222 % apply this to every coefficient in f. 223 % car cc is the degree list of this coefficient, 224 % cdr cc the coefficient itself. 225 % Finally collect the new pairs 226 % (degreelist . coefficient) 227 % 228 z := 229 for each cc in taycoefflist f join 230 for each cc2 in taylor2 (taycfsq cc, varlis, var0, n) 231 collect 232 taymakecoeff (append (taycfpl cc, taycfpl cc2), 233 taycfsq cc2); 234 % 235 % Append the new list to the Taylor template and return. 236 % 237 y := append(taytemplate f,list {varlis,var0,n,n+1}); 238 return make!-taylor!* (z, y, tayorig f, tayflags f) 239 end 240 % 241 % Last possible case: f is not a Taylor expression. 242 % Expand it. 243 % 244 else make!-taylor!* ( 245 taylor2 (f, varlis, var0, n), 246 list {varlis,var0,n,n+1}, 247 if !*taylorkeeporiginal then f else nil, 248 nil)$ 249 250symbolic procedure taylor2 (f, varlis, var0, n); 251 begin 252 scalar result, oldklist := get('taylor!*, 'klist); 253 result := errorset!* ({'taylor2e, mkquote f, mkquote varlis, mkquote var0, mkquote n}, nil); 254 resetklist('taylor!*, oldklist); 255 if atom result 256 then taylor!-error ('expansion, "(possible singularity!)") 257 else return car result 258 end$ 259 260symbolic procedure taylor2e (f, varlis, var0, n); 261 % 262 % taylor2e expands s.q. f w.r.t. varlis about var0 up to degree n. 263 % var is a list of kernels, which means that the expansion takes 264 % place in a homogeneous way if there is more than one kernel. 265 % The case that var0 is the special atom infinity has to be taken 266 % care of by the calling function(s). 267 % Expansion is carried out separately for numerator and 268 % denominator. This approach has the advantage of not producing 269 % complicated s.q.'s which usually appear if an s.q. is 270 % differentiated. The procedure is (roughly) as follows: 271 % if the denominator of f is free of var 272 % then expand the numerator and divide, 273 % else if the numerator is free of var expand the denominator, 274 % take the reciprocal of the Taylor series and multiply, 275 % else expand both and divide the two series. 276 % This fails if there are nontrivial dependencies, e.g., 277 % if a variable is declared to depend on a kernel in varlis. 278 % It is of course necessary that the denominator yields a unit 279 % in the ring of Taylor series. If not, an error will be 280 % signalled by invtaylor or quottaylor, resp. 281 % 282 if cdr varlis then taylor2hom (f, varlis, var0, n) 283 else if denr f = 1 then taylor2f (numr f, car varlis, var0, n, t) 284 else if not depends (denr f, car varlis) 285 then multintocoefflist (taylor2f (numr f, car varlis, var0, n, t), 286 1 ./ denr f) 287 else if numr f = 1 288 then delete!-superfluous!-coeffs 289 (invtaylor1 ({varlis,var0,n,n+1}, 290 taylor2f (denr f, car varlis, var0, n, nil)), 291 1, n) 292 else if not depends (numr f, car varlis) 293 then multintocoefflist 294 (invtaylor1 ({varlis,var0,n,n+1}, 295 taylor2f (denr f, car varlis, var0, n, nil)), 296 !*f2q numr f) 297 else begin scalar denom; integer n1; 298 denom := taylor2f (denr f, car varlis, var0, n, nil); 299 n1 := n + taymincoeff denom; 300 return 301 delete!-superfluous!-coeffs 302 (quottaylor1 ({varlis,var0,n1,n1+1}, 303 taylor2f (numr f, car varlis, var0, n1, t), 304 denom), 305 1, n) 306 end$ 307 308symbolic procedure taylor2f (f, var, var0, n, flg); 309 % 310 % taylor2f is the procedure that does the actual expansion 311 % of the s.f. f. 312 % It returns a list of the form 313 % ((deglis1 . coeff1) (deglis2 . coeff2) ... ) 314 % For the use of the parameter `flg' see below. 315 % 316 begin scalar x, y, z; integer k; 317 % 318 % Calculate once what is needed several times. 319 % var0 eq 'infinity is a special case that has already been taken 320 % care of by the calling functions by replacing var by 1/var. 321 % 322 if var0 eq 'infinity then var0 := 0; 323 x := list (var . var0); 324 y := simp list ('difference, var, var0); 325 % 326 % The following is a first attempt to treat expressions 327 % that have poles at the expansion point. 328 % Currently nothing more than factorizable poles, i.e. 329 % factors in the denominator are handled. 330 % We use the following trick to calculate enough terms: If the 331 % first l coefficients of the Taylor series are zero we replace n 332 % by n + 2l. This is necessary since we separately expand 333 % numerator and denominator of an expression. If the expansion of 334 % both parts starts with, say, the quadratic term we have to 335 % expand them up to order (n+2) to get the correct result up to 336 % order n. However, if the numerator starts with a constant term 337 % instead, we have to expand up to order (n+4). It is important, 338 % however, to drop the additional coefficients as soon as they are 339 % no longer needed. The parameter `flg' is used here to control 340 % this behaviour. The calling function must pass the value `t' if 341 % it wants to inhibit the calculation of additional coefficients. 342 % This is currently the case if the s.q. f has a denominator that 343 % may contain the expansion variable. Otherwise `flg' is used to 344 % remember if we already encountered a non-zero coefficient. 345 % 346 f := !*f2q f; 347 z := subs2 subsq (f, x); 348 if null numr z and not flg then n := n + 1 else flg := t; 349 y := list taymakecoeff ((list list 0), z); 350 k := 1; 351 while k <= n and not null numr f do 352 << f := multsq (diffsq (f, var), 1 ./ k); 353 % k is always > 0! 354 % subs2 added to simplify expressions involving roots 355 z := subs2 subsq (f, x); 356 if null numr z and not flg then n := n + 2 else flg := t; 357 if not null numr z then y := taymakecoeff(list list k, z) . y; 358 k := k + 1 >>; 359 return reversip y 360 end; 361 362symbolic procedure taylor2hom (f, varlis, var0, n); 363 % 364 % Homogeneous expansion of s.q. f wrt the variables in varlis, 365 % i.e. the sum of the degrees of the kernels is varlis is <= n 366 % 367 if null cdr varlis then taylor2e (f, list car varlis, var0, n) 368 else for each u in taylor2e (f, list car varlis, var0, n) join 369 for each v in taylor2hom (cdr u, 370 cdr varlis, 371 var0, 372 n - get!-degree taycfpl car u) 373 collect list (car taycfpl car u . taycfpl car v) . cdr v$ 374 375symbolic procedure taylorsamevar (u, varlis, var0, n); 376 begin scalar tp; integer mdeg, pos; 377 if cdr varlis 378 then taylor!-error ('not!-implemented, 379 "(homogeneous expansion in TAYLORSAMEVAR)"); 380 tp := taytemplate u; 381 pos := car var!-is!-nth (tp, car varlis); 382 tp := nth (tp, pos); 383 if taytpelpoint tp neq var0 384 then return taylor1 (if not null tayorig u then tayorig u 385 else simp!* preptaylor!* u, 386 varlis, var0, n); 387 mdeg := taytpelorder tp; 388 if n=mdeg then return u 389 else if n > mdeg 390 % 391 % further expansion required 392 % 393 then << lprim "Cannot expand further... truncated."; 394 return u >> ; 395 return make!-taylor!* ( 396 for each cc in taycoefflist u join 397 if nth (nth (taycfpl cc, pos), 1) > n 398 then nil 399 else list cc, 400 replace!-nth(taytemplate u,pos, 401 {varlis,taytpelpoint tp,n,n+1}), 402 tayorig u, tayflags u) 403 end$ 404 405 406COMMENT The `FULL' flag causes the whole term (including the 407 Taylor!* symbol) to be passed to SIMPTAYLOR!* ; 408 409symbolic procedure simptaylor!* u; 410 if taycoefflist u memq frlis!* or eqcar(taycoefflist u,'!~) 411 then !*tay2q u 412 else << 413 % 414 % Allow automatic combination of Taylor kernels. 415 % 416 if !*taylorautocombine and not taynomul!* and not ('taysimpsq memq mul!*) 417 then mul!* := aconc!* (mul!*, 'taysimpsq!-from!-mul); 418 !*tay2q resimptaylor u >>$ 419 420flag ('(taylor!*), 'full)$ 421 422put ('taylor!*, 'simpfn, 'simptaylor!*)$ 423 424COMMENT The following is necessary to properly process Taylor kernels 425 in substitutions; 426 427flag ('(taylor!*), 'simp0fn); 428 429endmodule; 430 431end; 432