Artifact 1a12be4c5927df2feb755df5fda81d3f145ebf80670dbb0dacb42e67fe4261d7:
- Executable file
r37/packages/plot/plotnum.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8883) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/plot/plotnum.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8883) [annotate] [blame] [check-ins using]
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;