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