Check-ins Using Download Hex

Artifact e7d428fc443adb136bf5a014e77b69328d291ee9f8907ca36c6f167796b86351:

module inttaylr;

% Author: James H. Davenport.

fluid '(const taylorasslist taylorvariable);

exports taylorform,taylorformp,taylorevaluate,return0,taylorplus,
         initialtaylorplus,taylorminus,initialtaylorminus,
         tayloroptminus,tayloroptplus,taylorctimes,initialtaylortimes,
         tayloroptctimes,taylorsqrtx,initialtaylorsqrtx,
         taylorquotient,initialtaylorquotient,taylorformersqrt,
         taylorbtimes,taylorformertimes,taylorformerexpt;

 symbolic procedure taylorform sq;
 if involvesf(denr sq,taylorvariable)
   then taylorformp list('quotient,tayprepf numr sq,tayprepf denr sq)
   else if 1 iequal denr sq
     then taylorformp tayprepf numr sq
     else taylorformp list('constanttimes,
                           tayprepf numr sq,
                           mk!*sq(1 ./ (denr sq)));
 % get division by a constant right.


 symbolic procedure taylorformp pf;
 if null pf
   then nil
   else if not dependsp(pf,taylorvariable)
     then taylorconst simp pf
     else begin
       scalar fn,initial,args,n;
       if atom pf
         then if pf eq taylorvariable
           then return taylorformp list ('expt,pf,1)
           else interr "False atom in taylorformp";
       % get 'x right as reduce shorthand for x**1.
       if taylorp pf
         then return pf;
       % cope with pre-expressed cases.
       % ***store-hack-1***
       % remove the (car pf eq 'sqrt) if more store is available.
       if (car pf eq 'sqrt) and
          (fn:=assoc(pf,taylorasslist))
         then go to lookupok;
       % look it up first.
       fn:=get(car pf,'taylorformer);
       if null fn
         then go to ordinary;
       fn:=apply1(fn,cdr pf);
       % ***store-hack-1***
       % remove the test if more store is available.
       if car pf eq 'sqrt
         then taylorasslist:=(pf.fn).taylorasslist;
       return fn;
       % cope with the special cases.
     ordinary:
       args := for each j in cdr pf collect taylorformp j;
       fn:=get(car pf,'tayloropt);
       if null fn
         then go to nooptimisation;
       fn:=apply1(fn,args);
       if fn
         then go to ananswer;
       % an optimisation has been made.
     nooptimisation:
       fn:=get(car pf,'taylorfunction);
       if null fn
         then interr "No Taylor function provided";
       fn:=fn.args;
       % fn is now the "how to compute" code.
       initial:=get(car pf,'initialtaylorfunction);
       if null initial
         then interr "No initial Taylor function";
       initial:=lispapply(initial,
                      list for each u in cdr fn collect firstterm u);
       % the first term in the expansion, or so we hope.
       n:=car initial;
       fn:=list(fn,n.n,initial);
       while null numr cdr initial do <<
             n:=n+1;
             if !*tra then lprim list("Increasing accuracy to",n);
             initial:=n.taylorevaluate(fn,n);
             fn:=list(car fn,n.n,initial);
             >>;
     ananswer:
       % ***store-hack-1***
       % uncomment this if more store is available;
       % taylorasslist:=(pf.fn).taylorasslist;
       return fn;
     lookupok:
       % These PRINT statements can be enabled in order to test the
       % efficacy of the association list
 %      printc "Taylor lookup succeeded";
 %      superprint car fn;
 %      printc length taylorasslist;
       return cdr fn
       end;


 symbolic procedure taylorevaluate(texpr,n);
 if n<taylorfirst texpr
   then nil ./ 1
   else if n>taylorlast texpr
     then tayloreval2(texpr,n)
     else begin
       scalar u;
       u:=assoc(n,taylorlist texpr);
       if u
         then return cdr u
         else return tayloreval2(texpr,n)
       end;


 symbolic procedure tayloreval2(texpr,n);
 begin
   scalar u;
   % actually evaluates from scratch.
   u:=apply3(taylorfunction texpr,n,texpr,cdr taylordefn texpr);
   if 'return0 eq taylorfunction texpr
     then return u;
   % no need to update with trivial zeroes.
   rplacd(cdr texpr,(n.u).taylorlist texpr);
   % update the association list.
   if n>taylorlast texpr
     then rplacd(taylornumbers texpr,n);
   % update the first/last pointer.
   return u
   end;


 symbolic procedure taylorconst sq;
 list('return0 . nil,0 . 0,0 . sq);


 symbolic procedure return0 (a,b,c);
 nil ./ 1;

 flag('(return0),'taylor);


 symbolic procedure firstterm texpr;
 begin
   scalar n,i;
   i:=taylorfirst texpr;
 trynext:
   n:=taylorevaluate(texpr,i);
   if numr n
     then return i.n;
   if i > 50
     then interr "Potentially zero Taylor series";
   i:=iadd1 i;
   rplaca(taylornumbers texpr,i);
   go to trynext
   end;


 symbolic procedure tayloroneterm u;
 % See if a Taylor expression has only one term.
  'return0 eq taylorfunction u and taylorfirst u=taylorlast u;


 % ***store-hack-1***;
 % uncomment this procedure if more store is available;
 % there is a smacro for this at the start of the file
 % for use if no store can be spared;
 %symbolic procedure tayshorten(save);
 %begin
 %  scalar z;
 %  % shortens the association list back to save,
 %    removing all the non-sqrts from it;
 %  while taylorasslist neq save do <<
 %    if caar taylorasslist eq 'sqrt
 %      then z:=(car taylorasslist).z;
 %    taylorasslist:=cdr taylorasslist >>;
 %  taylorasslist:=nconc(z,taylorasslist);
 %  return nil
 %  end;


 symbolic procedure tayprepf sf;
 if atom sf
   then sf
   else if atom mvar sf
     then taylorpoly makemainvar(sf,taylorvariable)
     else if null red sf
       then tayprept lt sf
       else list('plus,tayprept lt sf,tayprepf red sf);


 symbolic procedure tayprept term;
 if tdeg term = 1
   then if tc term = 1
     then tvar term
     else list('times,tvar term,tayprepf tc term)
   else if tc term = 1
     then list ('expt,tvar term,tdeg term)
     else list('times,list('expt,tvar term,tdeg term),
                    tayprepf tc term);


 symbolic procedure taylorpoly sf;
 % SF is a poly with MVAR = TAYLORVARIABLE.
 begin
   scalar tmax,tmin,u;
   tmax:=tmin:=ldeg sf;
   while sf do
     if atom sf or (mvar sf neq taylorvariable)
       then <<
         tmin:=0;
         u:=(0 . !*f2q sf).u;
         sf:=nil >>
       else <<
         u:=((tmin:=ldeg sf) . !*f2q lc sf) . u;
         sf:=red sf >>;
   return (list 'return0) . ((tmin.tmax).u)
   end;


 symbolic procedure taylorplus(n,texpr,args);
 mapply(function !*addsq,
        for each u in args collect taylorevaluate(u,n));


 symbolic procedure initialtaylorplus slist;
 begin
   scalar n,numlst;
   n:=mapply(function min2,mapovercar slist);
   % the least of the degrees.
   numlst:=nil;
   while slist do <<
     if caar slist iequal n
       then numlst:=(cdar slist).numlst;
     slist:=cdr slist >>;
   return n.mapply(function !*addsq,numlst)
   end;


 put ('plus,'taylorfunction,'taylorplus);
 put ('plus,'initialtaylorfunction,'initialtaylorplus);


 symbolic procedure taylorminus(n,texpr,args);
 negsq taylorevaluate(car args,n);


 symbolic procedure initialtaylorminus slist;
 (caar slist).(negsq cdar slist);


 put('minus,'taylorfunction,'taylorminus);
 put('minus,'initialtaylorfunction,'initialtaylorminus);


 flag('(taylorplus taylorminus),'taylor);


 symbolic procedure tayloroptminus(u);
 if 'return0 eq taylorfunction car u
   then taylormake(taylordefn car u,
                   taylornumbers car u,
                   taylorneglist taylorlist car u)
   else if 'taylorctimes eq taylorfunction car u
     then begin
       scalar const;
       u:=car u;
       const:=caddr taylordefn u;
       % the item to be negated.
       const:=taylormake(taylordefn const,
                         taylornumbers const,
                         taylorneglist taylorlist const);
       return taylormake(list(taylorfunction u,
                              argof taylordefn u,
                              const),
                         taylornumbers u,
                         taylorneglist taylorlist u)
       end
     else nil;
 put('minus,'tayloropt,'tayloroptminus);


 symbolic procedure taylorneglist u;
    for each v in u collect (car v . negsq cdr v);


 symbolic procedure tayloroptplus args;
 begin
   scalar ret,hard,u;
   u:=args;
   while u do <<
     if 'return0 eq taylorfunction car u
       then ret:=(car u).ret
       else hard:=(car u).hard;
     u:=cdr u >>;
   if null ret or
       null cdr ret
     then return nil;
   ret:=mapply(function joinret,ret);
   if null hard
     then return ret;
   rplaca(args,ret);
   rplacd(args,hard);
    return nil
   end;
 put('plus,'tayloropt,'tayloroptplus);


 symbolic procedure joinret(u,v);
 begin
   scalar nums,a,b,al;
   nums:=(min2(taylorfirst u,taylorfirst v).
          max2(taylorlast u,taylorlast v));
   al:=nil;
   u:=taylorlist u;
   v:=taylorlist v;
   for i:=(car nums) step 1 until (cdr nums) do <<
     a:=assoc(i,u);
     b:=assoc(i,v);
     if a
       then if b
         then al:=(i.!*addsq(cdr a,cdr b)).al
         else al:=a.al
       else if b
         then al:=b.al  >>;
   return taylormake(list 'return0,nums,al)
   end;




 % the operator constanttimes
 % has two arguments (actually a list)
 % 1) a form dependent on the taylorvariable
 % 2) a form which is not.


 % the operator binarytimes has two arguments (actually a list)
   % but behaves like times otherwise.


 symbolic procedure taylorctimes(n,texpr,args);
 !*multsq(taylorevaluate(car args,n-(taylorfirst cadr args)),
        taylorevaluate(cadr args,taylorfirst cadr args));


 symbolic procedure initialtaylortimes slist;
 % Multiply the variable by the constant.
 ((caar slist)+(caadr slist)). !*multsq(cdar slist,cdadr slist);


 symbolic procedure tayloroptctimes u;
 if 'taylorctimes eq taylorfunction car u
   then begin
     scalar reala,const,iconst,degg;
     % we have nested multiplication.
     reala:=argof taylordefn car u;
     % the thing to be multiplied by the two constants.
     const:=car taylorlist cadr u;
     %the actual outer constant: deg.sq.
     iconst:=caddr taylordefn car u;
     %the inner constant.
     degg:=(taylorfirst iconst)+(car const);
     iconst:=list(taylordefn iconst,
                   degg.degg,
                   degg.!*multsq(cdar taylorlist iconst,cdr const));
     return list('taylorctimes,reala,iconst).
                 ((((taylorfirst car u) + (car const)).
                         ((taylorlast car u) + (car const))).
		  for each j in taylorlist car u collect multconst j)
     end
   else if 'return0 eq taylorfunction car u
     then begin
       scalar const;
       const:=car taylorlist cadr u;
       % the actual constant:deg.sq.
       u:=car u;
       return (taylordefn u).
                   ((((taylorfirst u)+car const).
                         ((taylorlast u)+car const)).
		   for each j in taylorlist u collect multconst j)
       end
     else nil;


 symbolic procedure multconst v;
 % Multiplies v by const in deg.sq form.
 ((car v)+(car const)) . !*multsq(cdr v,cdr const);


 put('constanttimes,'tayloropt,'tayloroptctimes);
 put('constanttimes,'simpfn,'simptimes);
 put('constanttimes,'taylorfunction,'taylorctimes);
 put('constanttimes,'initialtaylorfunction,'initialtaylortimes);


 symbolic procedure taylorbtimes(n,texpr,args);
 begin
   scalar answer,n1,n2;
   answer:= nil ./ 1;
   n1:=car firstterm car args;
   % the first term in one argument.
   n2:=car firstterm cadr args;
   % the first term in the other.
   for i:=n1 step 1 until (n-n2) do
     answer:=!*addsq(answer,!*multsq(taylorevaluate(cadr args,n-i),
                                       taylorevaluate(car args,i)));
   return answer
   end;




 put('binarytimes,'taylorfunction,'taylorbtimes);
 put('binarytimes,'initialtaylorfunction,'initialtaylortimes);
 put('binarytimes,'simpfn,'simptimes);


symbolic procedure taylorformertimes arglist;
begin
  scalar const,var,degg,wsqrt,negcount;  % u;
  negcount:=0;
  degg:=0;% the deggrees of any solitary x we may meet.
  const:=nil;
  var:=nil;
  wsqrt:=nil;
  while arglist do <<
    if dependsp(car arglist,taylorvariable)
      then if and(eqcar(car arglist,'expt),
                        cadar arglist eq taylorvariable,
                        numberp caddar arglist)
        then degg:=degg+caddar arglist
% removed JHD 21.8.86 - while it is anoptimisation,
% it runs the risk of proving that -1 = +1 by ignoring the
% number of "i" needed - despite the attempts we went to.
%        else if eqcar(car arglist,'sqrt)
%          then <<
%            u:=argof car arglist;
%            wsqrt := u . wsqrt;
%            if minusq cdr firstterm taylorformp u
%              then negcount:=1+negcount >>
          else if car arglist eq taylorvariable
            then degg:=degg + 1
            else var:=(car arglist).var
      else const:=(car arglist).const;
    arglist:=cdr arglist >>;
  if wsqrt
    then if cdr wsqrt
      then var:=list('sqrt,prepsq simptimes wsqrt).var
      else var:=('sqrt.wsqrt).var;
  if var
    then var:=mapply(function (lambda u,v;
                               list('binarytimes,u,v)),var);
  % insert binary multiplications.
  negcount:=negcount/2;
  if onep cdr divide(negcount,2)
    then const:= (-1).const;
  % we had an odd number of (-1) from i*i.
  if const or (degg neq 0)
    then <<
      if const
        then const:=simptimes const
        else const:=1 ./ 1;
      const:=taylormake(list 'return0,degg.degg,list(degg.const));
      if null var
        then var:=const
        else var:=list('constanttimes,var,const) >>;
  return taylorformp var
  end;

put('times,'taylorformer,'taylorformertimes);




flag('(taylorbtimes taylorctimes taylorquotient),'taylor);
symbolic procedure taylorformerexpt arglist;
begin
  scalar base,expon;
  base:=car arglist;
  expon:=simpcar cdr arglist;
  if (denr expon neq 1) or
     (not numberp numr expon)
    then interr "Hard exponent";
  expon:=numr expon;
  if base neq taylorvariable
    then interr "Hard base";
  return list('return0 . nil,expon.expon,expon.(1 ./ 1))
  end;
put ('expt,'taylorformer,'taylorformerexpt);


symbolic procedure initialtaylorquotient slist;
(caar slist - caadr slist). !*multsq(cdar slist,!*invsq cdadr slist);


symbolic procedure taylorquotient(n,texpr,args);
begin
  % problem is texpr=b/c or c*texpr=b.
  scalar sofar,b,c,cfirst;
  b:=car args;
  c:=cadr args;
  cfirst:=taylorfirst c;
  sofar:=taylorevaluate(b,n+cfirst);
  for i:=taylorfirst texpr step 1 until n-1 do
    sofar:=!*addsq(sofar,!*multsq(taylorevaluate(texpr,i),
                              negsq taylorevaluate(c,n+cfirst-i)));
  return !*multsq(sofar,!*invsq taylorevaluate(c,cfirst))
  end;


put('quotient,'taylorfunction,'taylorquotient);
put('quotient,'initialtaylorfunction,'initialtaylorquotient);


% symbolic procedure minusq sq;
%    if null sq then nil
%     else if minusf numr sq then not minusf denr sq
%     else minusf denr sq;

% This is wrapped round TAYLORFORMERSQRT2 in order to
% remove the innards of the SQRT from the asslist.
% note the precautions for nested SQRTs.

symbolic procedure taylorformersqrt arglist;
% ***store-hack-1***;
% Uncomment these lines if more store is available.
%begin
%  scalar z;
%  z:=taylorasslist;
%  if sqrtsintree(car arglist,taylorvariable)
%    then return taylorformersqrt2 arglist;
%  arglist:=taylorformersqrt2 arglist;
%  taylorasslist:=z;
%  return arglist
%  end;
%
%
%symbolic procedure taylorformersqrt2 arglist;
begin
  scalar f,realargs,ff,realsqrt;
  realargs:=taylorformp carx(arglist,'taylorformersqrt2);
  f:=firstterm realargs;
  if not evenp car f
    then interr "Extra sqrt substitution needed";
  if and(0 iequal car f,
         1 iequal numr cdr f,
         1 iequal denr cdr f)
    then return taylorformp list('sqrtx,realargs);
  % if it starts with 1 already then it is easy.
  ff:=- car f;
  ff:=list(list 'return0,ff.ff,ff.(!*invsq cdr f));
  % ff is the leading term in the expansion of realargs.
  realsqrt:=list('sqrtx,list('constanttimes,realargs,ff));
  ff:=(car f)/2;
  return taylorformp list('constanttimes,
                          realsqrt,
                          list(list 'return0,
                               ff.ff,
                               ff.(simpsqrtsq cdr f)))
  end;


put('sqrt,'taylorformer,'taylorformersqrt);


symbolic procedure initialtaylorsqrtx slist;
0 . (1 ./ 1);
% sqrt(1+ ...) = 1+....


symbolic procedure taylorsqrtx(n,texpr,args);
begin scalar sofar;
  sofar:=taylorevaluate(car args,n);
  % (1+.....+a(n)*x**n)**2
  % = ....+x**n*(2*a(n)+sum(0<i<n,a(i)*a(n-i))).
  % So a(n)=(coeff(x**n)-sum) /2.
  for i:=1 step 1 until (n-1) do
    sofar:=!*addsq(sofar,negsq !*multsq(taylorevaluate(texpr,i),
                                    taylorevaluate(texpr,n-i)));
  return multsq(sofar,1 ./ 2)
  end;


flag('(taylorsqrtx),'taylor);
put('sqrtx,'taylorfunction,'taylorsqrtx);
put('sqrtx,'initialtaylorfunction,'initialtaylorsqrtx);

endmodule;

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]