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