reduce-historical/r38/packages/plot/plotnum.red

module plotnum; % Numeric evaluation of algebraic expressions.

fluid '(plotsynerr!* ploteval!-alist2!*);

global '(!*plotinterrupts);

flag('(plus plus2 difference times times2 quotient exp expt expt!-int
%      minus sin cos tan cot asin acos acot atan log),'plotevallisp);
       minus sin cos tan cot asin acos acot atan abs log),'plotevallisp);

symbolic procedure plotrounded d;
 % Switching rounded mode safely on and off.
  begin scalar o,!*protfg,c,r,!*msg;
   !*protfg:=t;
   if null d then
    <<c:=!*complex; r:=!*rounded;
      if c then <<setdmode('complex,nil); !*complex := nil>>;
      if not r and dmode!* then
          <<o:=get(dmode!*,'dname); setdmode(o,nil)>>;
      o:={o,r,!*roundbf,c,precision 0};
      if not r then <<!*rounded:=t; setdmode('rounded,t)>>;
      !*roundbf:=nil;
      if c then <<setdmode('complex,t); !*complex := t>>;
      return o;
    >>
   else
    <<
       % reconstruct the previous configuration.
      if !*complex then setdmode('complex,nil);
      if null(!*rounded:=cadr d) then setdmode('rounded,nil);
      !*roundbf:=caddr d; 
      if car(d) then setdmode(car d,t);
      if !*complex then setdmode('complex,t);
      precision car cddddr d;
    >>;
   end;

symbolic procedure adomainp u;
 % numberp test in an algebraic form.
   numberp u or (pairp u and idp car u and get(car u,'dname))
             or eqcar(u,'minus) and adomainp cadr u;

symbolic procedure revalnuminterval(u,num);
 % Evaluate u as interval; numeric bounds required if num=T.
  begin scalar l;
    if not eqcar(u,'!*interval!*) then typerr(u,"interval");
    l:={reval cadr u,reval caddr u};
    if null num or(adomainp car l and adomainp cadr l)then return l;
    typerr(u,"numeric interval");
  end;

ploteval!-alist2!*:={{'x . 1},{'y . 2}};

symbolic procedure plot!-form!-prep(f,x,y);
  % f is a REDUCE function expression depending of x and y. 
  % Compute a form which allows me to evaluate "f(x,y)" as
  % a LISP expr.
  {'lambda,'(!&1 !&2),
     {'plot!-form!-call2,mkquote rdwrap f,mkquote f,
                mkquote x,'!&1,
                mkquote y,'!&2}};

symbolic procedure plot!-form!-call2(ff,f,x,x0,y,y0);
 % Here I hack into the statically allocated a-list in order
 % to save cons cells.
  begin scalar a;
     a:=car ploteval!-alist2!*;
     car a := x; cdr a := x0;
     a:=cadr ploteval!-alist2!*;   
     car a:= y; cdr a:= y0;
     return plotevalform(ff,f,ploteval!-alist2!*);
  end;

% The following routines are used to transform a REDUCE algebraic
% expression into a form which can be evaluated directly in LISP.

symbolic procedure rdwrap f;
  begin scalar r,!*msg,!*protfg;
    !*protfg:=t;
    r:=errorset({'rdwrap1,mkquote f},nil,nil);
    return if errorp r then 'failed else car r;
  end;

symbolic procedure rdwrap1 f;
if numberp f then float f
  else if f='pi then 3.141592653589793238462643
  else if f='e then 2.7182818284590452353602987
  else if f='i and !*complex then rederr "i not LISP evaluable" 
  else if atom f then f
  else if get(car f,'dname) then rdwrap!-dm f
  else if eqcar(f, 'MINUS) then
    begin scalar x;
       x := rdwrap1 cadr f;
       return if numberp x then minus float x else {'MINUS, x}
    end
  else if eqcar(f,'expt) then rdwrap!-expt f
  else
   begin scalar a,w;
     if null getd car f or not flagp(car f,'plotevallisp)
                then typerr(car f,"Lisp arithmetic function (warning)");
     a:=for each c in cdr f collect rdwrap1 c;
     if (w:=atsoc(car f,'((plus.plus2)(times.times2)))) 
             then return rdwrapn(cdr w,a);
     return car f . a;
   end;

symbolic procedure rdwrapn(f,a);
  % Unfold a n-ary arithmetic call.
   if null cdr a then car a else {f,car a,rdwrapn(f,cdr a)};

symbolic procedure rdwrap!-dm f;
 % f is a domain element.
  if car f eq '!:RD!: then 
          if atom cdr f then cdr f else bf2flr f
  else if car f eq '!:CR!: then rdwrap!-cr f
  else if car f eq '!:GI!:
   then rdwrap!-cmplex(f,float cadr f,float cddr f)
  else if car f eq '!:DN!: then rdwrap2 cdr f
  else << plotsynerr!*:=t;
       rerror(PLOTPACKAGE, 32, {get(car f, 'DNAME),
                                "illegal domain for PLOT"})
    >>;
  
symbolic procedure rdwrap!-cr f;
  begin scalar re,im;
    re := cadr f;
    if not atom re then re := bf2flr re;
    im := cddr f;
    if not atom im then im := bf2flr im;
    return rdwrap!-cmplx(f,re,im);
  end;

symbolic procedure rdwrap!-cmplx(f,re,im);
   if abs im * 1000.0 > abs re then typerr(f,"real number") else re;

symbolic procedure plotrepart u;
   if car u eq '!:rd!: then u else
   if car u memq '(!:gi!: !:cr!:) then '!:rd!: . cadr u;

symbolic procedure rdwrap!-expt f;
  % preserve integer second argument.
  if fixp caddr f then 
     if caddr f>0 then {'expt!-int,rdwrap1 cadr f,caddr f}
       else {'quotient,1,{'expt!-int,rdwrap1 cadr f,-caddr f}}
    else {'expt,rdwrap1 cadr f, rdwrap caddr f};

symbolic procedure rdwrap2 f;
% Convert from domain to LISP evaluable value.
  if atom f then f else float car f * 10^cdr f;

symbolic procedure rdwrap!* f;
  % convert a domain element to float.
  if null f then 0.0 else rdwrap f;

symbolic procedure rdunwrap f; 
    if f=0.0 then 0 else if f=1.0 then 1 else '!:rd!: . f;

symbolic procedure expt!-int(a,b); expt(a,fix b);

symbolic procedure plotevalform(ff,f,a);
  % ff is LISP evaluable,f REDUCE evaluable.
  begin scalar u,w,!*protfg,mn,r,!*msg;
    !*protfg := t;
%   if ff='failed then goto non_lisp;
    if ff eq 'failed or not atom ff and 'failed memq ff
      then goto non_lisp;
% WN 12. May.97  plot(e^(abs x)); bombed
    u:= errorset({'plotevalform1,mkquote ff,mkquote a},nil,nil);
    if not ploterrorp u and numberp (u:=car u) then 
    <<if abs u > plotmax!* then return plotoverflow plotmax!* else
      return u;
    >>;
 non_lisp:
    w := {'simp, mkquote
     sublis(
       for each p in a collect (car p.rdunwrap cdr p),
       f)};
    u := errorset(w,nil,nil);
    if ploterrorp u or (u:=car u) eq 'failed then return nil;
    if denr u neq 1 then return nil;
    u:=numr u;
    if null u then return 0; % Wn
    if numberp u then <<r:=float u; goto x>>;
    if not domainp u or not memq(car u,'(!:rd!: !:gi!: !:cr!:)) 
          then return nil;
    if evalgreaterp(plotrepart u, fl2rd plotmax!*) then
     return plotoverflow plotmax!* else
    if evalgreaterp(fl2rd (-plotmax!*),plotrepart u) then
     return plotoverflow (-plotmax!*);
    r:=errorset({'rdwrap!-dm,mkquote u},nil,nil);
    if errorp r or(r:=car r) eq 'failed then return nil;
  x: return if mn then -r else r;
  end;

symbolic procedure plotoverflow u;
   <<if not !*plotoverflow then 
      lprim "plot number range exceeded";
     !*plotoverflow := t;
     'overflow . u
   >> where !*protfg = nil;

symbolic procedure plotevalform0(f,a);
  (if ploterrorp u 
       then typerr(f,"expression for plot") else car u)
     where u=
   errorset({'plotevalform1,mkquote f,mkquote a},nil,nil);

%symbolic procedure plotevalform1(f,a);
% begin scalar x,w;
%  if numberp f then return float f;
%  if (x:=assoc(f,a)) then return plotevalform1(cdr x,a);
%  if not atom f and flagp(car f,'plotevallisp) then
%    return eval 
%      (car f . for each y in cdr f collect plotevalform1(y,a));
%  if not atom f then f :=
%    car f . for each y in cdr f collect prepf plotevalform2(y,a);
%  w:=simp f;
%  if denr w neq 1 or not domainp numr w then goto err;
%  return rdwrap!* numr w;
% err:
%  plotsynerr!*:=t;
%  typerr(prepsq w,"numeric value");
%end;

symbolic procedure plotevalform1(f,a);
 begin scalar x;
  if numberp f then return float f;
  if (x:=assoc(f,a)) then return plotevalform1(cdr x,a);
  if atom f then go to err;
  return if cddr f then 
     idapply(car f,{plotevalform1(cadr f,a),plotevalform1(caddr f,a)})
    else
     idapply(car f,{plotevalform1(cadr f,a)}); 
 err:
  plotsynerr!*:=t;
  typerr(prepsq f,"numeric value");
end;


%symbolic procedure plotevalform2(f,a);
% begin scalar x,w;
%  if fixp f then return f;
%  if floatp f then return rdunwrap f;
%  if (x:=assoc(f,a)) then return plotevalform2(cdr x,a);
%  if not atom f and flagp(car f,'plotevallisp) then
%    return rdunwrap eval 
%      (car f . for each y in cdr f collect plotevalform1(y,a));
%  if not atom f then f :=
%    car f . for each y in cdr f collect prepf plotevalform2(y,a);
%  w:=simp f;
%  if denr w neq 1 or not domainp numr w then goto err;
%  return numr w;
% err:
%  plotsynerr!*:=t;
%  typerr(prepsq w,"numeric value");
%end;

symbolic procedure ploterrorp u;
   if u member !*plotinterrupts 
      then rederr prin2 "***** plot killed"
    else atom u or cdr u;

endmodule;

end;