reduce-historical/r38/packages/numeric/numeval.red

module numeval; % numeric evaluation of algebraic expressions.

% Control of accuracy and precision.
%
%  precision:    number of digits used for computation,
%  accuracy:     target precision of the results; 
%
%  precision might be modified automatically for reaching the
%  desired accuracy.

symbolic procedure accuracycontrol(u,da,di);
  % u is an evaluated parameter list. Extract
  % accuracy and iterations. If not present, take
  % given default values.
  begin scalar x,n,v;
   v:=u;
   accuracy!*:=da; iterations!*:=di;
   while v do
   <<x:=car v; v:= cdr v;
     if eqcar(x,'equal) and cadr x memq'(accuracy iterations) then
     <<u:=delete(x,u); n:=caddr x;
       if cadr x='accuracy then accuracy!*:=n
           else iterations!*:=n;
   >>>>;
   return u;
 end;

symbolic procedure update!-precision l;
  % l is a list of domain elements. IF the minimum of their absolute
  % values is smaller than 10**(precision*+3), the precision is
  % increased.
  begin scalar mn; integer zp,p;
   mn:=update!-precision1 l; % positive minimum.
   if null mn then return nil;
   p := precision 0; zp:=!:recip expt(10,p-3);
   dm!: <<
     if mn>zp then return nil;
     while mn<zp do <<p:=p+1;zp:=zp/10>>;
     >>;
   precmsg p;
  end;

symbolic procedure update!-precision1 l;
 dm!:( begin scalar x,y,z;
      while l do
      <<x:=car l; l:=cdr l;
        if not zerop x then
        <<y:=abs x; z:=if null z or y<z then y else z>>
      >>;
      return z;
     end );

% Switching of mode and introduction of specific simplifiction
% rules.

algebraic<<

  rules_rd :=
   {~u ** ~x => exp(log u * x) when not numberp x};

  procedure switch!-mode!-rd!-alg u;
   if u=0 then
     <<for all x clear exp x;
       let rules_rd;
     >> else <<
       for all x let exp x=e**x;
       clearrules rules_rd; >>;

>>;

symbolic procedure switch!-mode!-rd u;
  begin scalar oldmode,prec,ne;
  if null u then
    <<if not memq(dmode!*,'(!:rd!: !:cr))then
       <<oldmode:=t; setdmode('rounded,t)>>;
     ne := !*noequiv;
     !*noequiv:=t; prec:=precision 0;
     switch!-mode!-rd!-alg 0;
     return list(oldmode,prec,!*roundbf,ne)
    >> else <<
     if car u then setdmode('rounded,nil);
     precision cadr u;
     !*roundbf := caddr u;
     !*noequiv := cadddr u;
     switch!-mode!-rd!-alg 1;
    >>;
  end;

% Support for the numerical (=domain specific) evaluation of 
% algebraic equations.
%
% Roughly equivalent:
%   evaluate(u,v,p) = numr subsq(simp u,pair(v,p))
% but avoiding resimplification as we know that all components
% must evaluate to domain elements.

fluid '(dmarith!* !*evaluateerror);

dmarith!*:= '(
   (difference . !:difference) (quotient . !:!:quotient)
   (minus . negf) (sqrt . num!-sqrtf)
   (expt . !:dmexpt) (min . dm!:min) (max . dm!:max));

symbolic procedure evaluate(u,v,p);
   begin scalar a,r,!*evaluateerror,msg;
    msg := not !*protfg;
    a:= pair(v,p);
    r := errorset(list('evaluate0,mkquote u,mkquote a),msg,nil)
           where !*msg=nil,!*protfg=t;
    if errorp r then rederr
        "error during function evaluation (e.g. singularity)";
    return car r;
  end;

symbolic procedure evaluate!*(u,v,p);
  % same as evaluate, but returning arithmetic (but not syntactical)
  % errors to the caller like errorset.
   begin scalar a,r,!*evaluateerror;
    a:= pair(v,p);
    r := errorset(list('evaluate0,mkquote u,mkquote a),nil,nil)
          where !*msg=nil,!*protfg=t;
    erfg!* := nil;
    if null !*evaluateerror then return r else
         evaluate0(u,a); % once more, but unprotected.
  end;

symbolic procedure evaluate0(u,v);
   evaluate1(evaluate!-horner u,v);

symbolic procedure evaluate1(u,v);
    if numberp u or null u then force!-to!-dm u else
    if pairp u and get(car u,'dname) then u else
     (if a then cdr a else
      if atom u then
             if u='i then
              if (u:=get(dmode!*,'ivalue)) then numr apply(u,'(nil))
                else rederr "i used as indeterminate value"
             else if u='e or u='pi then force!-to!-dm numr simp u
             else <<!*evaluateerror:=t; typerr(u,"number")>>
             else evaluate2(car u,cdr u,v)
                ) where a=assoc(u,v);

symbolic procedure evaluate2(op,pars,v);
 if op='!:dn!: then numr dn!:simp pars else 
   <<pars:=for each p in pars collect evaluate1(p,v);
          % arithmetic operator.
     if op='plus then !:dmpluslst pars else
     if op='times then !:dmtimeslst pars else
     if(a:=assoc(op,dmarith!*)) then apply(cdr a,pars) else
          % elementary function, applicable for current dmode.
     if pairp car pars and (a:=get(op,caar pars)) then
          apply(a,pars) else
	 % apply REDUCE simplifier otherwise
     force!-to!-dm numr simp (op . for each p in pars collect prepf p)
    >> where a=nil;

symbolic procedure list!-evaluate(u,v,p);
     for each r in u collect evaluate(r,v,p);

symbolic procedure matrix!-evaluate(u,v,p);
     for each r in u collect list!-evaluate(r,v,p);

symbolic procedure !:dmexpt(u,v);
    (if fixp n then !:expt(u,n) else
      <<u:=force!-to!-dm u; v:=force!-to!-dm v;
        if !*complex then crexpt!*(u,v) else rdexpt!*(u,v)
      >>) where n=!:dm2fix v;

symbolic procedure dm!:min(a,b);
   dm!:(if a > b then b else a);

symbolic procedure dm!:max(a,b);
   dm!:(if a > b then a else b);

symbolic procedure !:dm2fix u;
  if fixp u then u else
    (if fixp(u:=int!-equiv!-chk u) then u else
     if null u then 0 else
     if floatp cdr u and 0.0=cdr u-fix cdr u then fix cdr u
     else u) where !*noequiv=nil;

symbolic procedure !:dmtimeslst u;
    if null u then 1 else
    if null cdr u then car u else
       dm!:(car u * !:dmtimeslst cdr u);

symbolic procedure !:dmpluslst u;
    if null u then 0 else
    if null cdr u then car u else
       dm!:(car u + !:dmpluslst cdr u);

symbolic procedure !:!:quotient(u,v);
   !:quotient(u,if fixp v then i2rd!* v else v);

% vector/matrix arithmetic

symbolic procedure list!+list(l1,l2);
   if null l1 then nil else
       dm!: (car l1 + car l2) . list!+list(cdr l1,cdr l2);

symbolic procedure list!-list(l1,l2);
   if null l1 then nil else
       dm!: (car l1 - car l2) . list!-list(cdr l1,cdr l2);

symbolic procedure scal!*list(s,l);
  if null l then nil else
     dm!: (s * car l) . scal!*list(s,cdr l) ;

symbolic procedure innerprod(u,v);
     if null u then 0 else
     dm!: ( car u * car v + innerprod(cdr u,cdr v) );

symbolic procedure conjlist u;
dm!:(if not !*complex then u else
    for each x in u collect
       repartf x - numr apply(get(dmode!*,'ivalue),'(nil))*impartf x );

symbolic procedure normlist u;
  dm!:(sqrt innerprod(u, conjlist u));

symbolic procedure mat!*list(m,v);
    if null cdr m then scal!*list(car v,car m) else
    for each r in m collect innerprod(r,v);

symbolic procedure num!-sqrtf a;
   if domainp a then
       apply(get('sqrt,dmode!*),list force!-to!-dm a)
     else
     <<a:=simpsqrt list prepf a;
       if not onep denr a or not domainp numr a
         then rederr "sqrtf called in wrong mode"
       else numr a>>;

symbolic procedure force!-to!-dm a;
   if not domainp a then typerr(prepf a,"number") else
   if null a then force!-to!-dm 0 else
   if numberp a then apply(get(dmode!*,'i2d),list a) else
   if pairp a and car a = dmode!* then a else
    (if fcn then apply(fcn,list a) else
       rederr list("conversion error with ",a)
     ) where fcn=(pairp a and get(car a,dmode!*));

symbolic procedure printsflist(x);
   begin integer n;
    writepri("(",nil);
    for each y in x do
       <<if n>0 then writepri(" , ",nil);
         n:=n+1;
         writepri(mkquote prepf y,nil)>>;
    writepri(")",nil);
  end;

fluid '(horner!-cache!*);

symbolic procedure evaluate!-horner u;
  (if w then cdr w else
   <<w:=simp!* u; 
     w:=prepsq(hornerf numr w ./ hornerf denr w);
     horner!-cache!* := (u.w). horner!-cache!*;
     w
   >>) where w=assoc(u,horner!-cache!*);

endmodule;

end;