module numeric;  % Header module for the numeric package and
                 % support of numerical evaluation of symbolic
                 % expressions.

% Author: Herbert Melenk.

% Copyright (c) 1993 ZIB Berlin, RAND.  All rights reserved.

create!-package('(numeric numeval numsolve newton steepstd bounds
                  numint numfit chebysh rungeku),
                '(contrib numeric));

fluid '(!*noequiv);
fluid '(accuracy!*);
global '(iterations!* !*trnumeric);
switch trnumeric;

% Create .. as the infix operator if not yet done.

newtok '( (!. !.) !*interval!*);
precedence .., or;
algebraic operator ..;
put('!*interval!*,'prtch,'! !.!.! );

% some common utilities

% intervals

symbolic procedure adomainp u;
   numberp u or (pairp u and idp car u and get(car u,'dname));

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;

symbolic procedure mkinterval(u,v);
   list('!*interval!*,u,v);

% Easy coding of numerical procedures with REDUCE:
%
%  In statements or procedure bodies tagged with "dm:" all 
%  arithmetic function calls are replaced by REDUCE domain 
%  arithmetic.

symbolic macro procedure dm!: u;
  subla('((plus2 . !:plus)(times2 . !:times)
          (plus . !:plusn)(times . !:timesn)
          (quotient . !:!:quotient)
          (difference . !:difference)
          (minus . !:minus)
          (minusp . !:minusp)
          (zerop . !:zerop)
          (lessp . (lambda(a b)(!:minusp (!:difference a b))))
          (greaterp . (lambda(a b)(!:minusp (!:difference b a))))
          (leq . (lambda(a b)(not (!:minusp (!:difference b a)))))
          (geq . (lambda(a b)(not (!:minusp (!:difference a b)))))
          (sqrt . num!-sqrtf)
          (abs . absf)
          (min2 . dm!:min)
          (max2 . dm!:max)
          (min . dm!:min)
          (max . dm!:max)
         ) , cadr u);

%wrappers for n-ary plus and times

symbolic macro procedure !:plusn u;
  if null cddr u then cadr u else
   list('!:plus,cadr u,'!:plusn . cddr u);

symbolic macro procedure !:timesn u;
  if null cddr u then cadr u else
   list('!:times,cadr u,'!:timesn . cddr u);

endmodule;


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;
  if null u then
    <<if not memq(dmode!*,'(!:rd!: !:cr))then
       <<oldmode:=t; setdmode('rounded,t)>>;
     !*noequiv:=t; prec:=precision 0;
     switch!-mode!-rd!-alg 0;
     return list(oldmode,prec,!*roundbf)
    >> else <<
     if car u then setdmode('rounded,nil);
     precision cadr u;
     !*roundbf := caddr 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;
    a:= pair(v,p);
    r := errorset(list('evaluate1,mkquote u,mkquote a),t,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('evaluate1,mkquote u,mkquote a),nil,nil)
          where !*msg=nil,!*protfg=t;
    erfg!* := nil;
    if null !*evaluateerror then return r else
         evaluate1(u,a); % once more, but unprotected.
  end;

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);
   <<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 siplificator 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;

endmodule;


module numsolve;  % root finding.

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1991, all rights resrved
% Version 2: termination and damping criterion modified to affine
%            invariance.

fluid '(!*noequiv accuracy!* !*exptexpand);
global '(iterations!* !*trnumeric erfg!*);

symbolic procedure rdsolveeval u;
     % interface function;
  begin scalar e,l,vars,y,z,p,r,erfg,mode,all,!*exptexpand;
        integer n;
    erfg:= erfg!*; erfg!* := nil;
    u := for each x in u collect reval x;
    u:=accuracycontrol(u,6,50);
    if (all:='all memq u) then u:=delete('all,u);
    e :=car u; u :=cdr u;
      % construct the function vector.
    e:=if eqcar(e,'list) then cdr e else list e;
    e := for each f in e collect
     if eqexpr (f:=reval f) then !*eqn2a f else f;
    n := length e;
      % starting point & variables.
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if eqcar(x,'equal) then
       <<z:=cadr x; y:=caddr x;
         if eqcar(y,'!*interval!*) then mode:='i;
       >> else <<z:=x; y:=random 100>>;
       vars:=z . vars; p := y . p>>;
    vars := reversip vars; p := reversip p;
    if not(n=length vars) then
    <<
   %  lprim "numbers of variables and functions don't match;";
   %  lprim "using steepest descent method instead of Newton";
   %  return rdsolvestdeval list ('list.e,'list.u,
   %                      {'equal,'accuracy,accuracy!*},
   %                      {'equal,'iterations,iterations!*});
      rederr "numbers of variables and functions don't match"
    >>;
    if mode='i and length vars>1 or mode neq 'i and all then
           rederr "mode for num_solve not implemented";
    r:=if mode='i then rd!-solve!-interv(e,vars,p,n,all)
        else rdnewton0(e,vars,p,n);
    erfg!* := erfg;
    return r;
  end;

put('num_solve,'psopfn,'rdsolveeval);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  finding root in an interval (secant method and related)
%

symbolic procedure rd!-solve!-interv(e,vars,p,n,all);
 (begin scalar u,fu,l,fl,x,acc,r,de,w,oldmode;
    if length p>1 then typerr('list . p,"single interval");
    p:=car p; e:=car e; x:=car vars;
    w := {'boundseval,mkquote {e,{'equal,x,p}}};
   if not memq(dmode!*,'(!:rd!: !:cr!:))then
       <<oldmode:=t; setdmode('rounded,t)>>;
    if errorp errorset(w,nil,nil) then
     typerr(e,"bounded expression");
    acc := !:!:quotient(1,expt(10,accuracy!*));
    l:=evaluate(cadr p,nil,nil); u:=evaluate(caddr p,nil,nil);
    fl:=evaluateuni(e,x,l); fu:=evaluateuni(e,x,u);
    if not all then
    dm!:(if zerop fl then return prepf l else
         if zerop fu then return prepf u else
         if fl*fu<0 then
          r := regula!-falsi(e,x,l,fl,u,fu,acc));
    if null r then de := reval {'df,e,x};
    if null r and not all then
       r:=rd!-solve!-trynewton(e,de,x,l,fl,u,fu,acc);
    if null r then
       r:=rd!-solve!-hardcase(e,de,x,l,fl,u,fu,acc,all);
    if oldmode then setdmode('rounded,nil);
    if r then return r;
    rederr "not applicable";
   end) where !*roundbf=!*roundbf;

symbolic procedure make!-rdsolve!-result1(x,u);
   'list.for each r in u collect {'equal,x,prepf r};

symbolic procedure rd!-solve!-trynewton(e,de,x,l,fl,u,fu,acc);
  begin scalar r,invde;
    invde := {'quotient,1,de};
  dm!: <<
    if fu*evaluateuni(de,x,u) >0 then
      r := rdnewton2({e},{{invde}},{x},acc,{u},'root,l,u);
    if null r and fl*evaluateuni(de,x,l) <0 then
      r := rdnewton2({e},{{invde}},{x},acc,{l},'root,l,u);
    if r and (r:=car r) and l <= r and r <= u then
         r := make!-rdsolve!-result1(x,{r}) else r:=nil;
     >>;
    return r;
  end;

symbolic procedure regula!-falsi(e,x,l,fl,u,fu,acc);
   make!-rdsolve!-result1(x,
       {regula!-falsi1(e,x,l,fl,u,fu,acc,0,1)});

symbolic procedure regula!-falsi1(e,x,l,fl,u,fu,acc,mode1,mode2);
  % modified regula falsi: using bisection in order to
  % traverse root as often as possible.
 dm!: begin scalar m,fm;
     if mode1=mode2 then m:=(u+l)/2
        else m := l*fu/(fu-fl) + u*fl/(fl-fu);
     if (u-l) < acc then return m;
     fm := evaluateuni(e,x,m);
     if zerop fm then return m;
     return if fl*fm<0
       then regula!-falsi1(e,x,l,fl,m,fm,acc,nil,mode1)
       else regula!-falsi1(e,x,m,fm,u,fu,acc,t,mode1);
    end;

symbolic procedure mkminus u; {'minus,u};

symbolic procedure evaluateuni(e,x,p);
     evaluate(e,{x},{p});

symbolic procedure dmboundsuni(e,x,l,u);
  begin scalar i;
     i:=boundseval {e,{'equal,x,{'!*interval!*,prepf l,prepf u} }};
     return numr simp cadr i . numr simp caddr i;
  end;

symbolic procedure rd!-solve!-hardcase(e,de,x,l,fl,u,fu,acc,all);
 dm!: begin scalar il1,il2,pt,pt1,iv,r,rr,b,b1,q,d;
     integer lev;
   il1:={(l.fl) .(u.fu)};
  main_loop:
    lev:=lev+1; il2:=nil;
    if null il1 then goto done;
    il1 := reverse il1;
    rd!-solve!-hardcase!-print(lev,il1);
  loop:
    if null il1 then goto bottom;
    iv :=car il1; il1:= cdr il1;
    l:=caar iv; fl:=cdar iv; u:=cadr iv; fu:=cddr iv;
    if (d:=u-l) <0 then goto loop; %empty
    if abs fl<acc then<<pt:=l;goto root_found>>;
    if abs fu<acc then<<pt:=u;goto root_found>>;
    b:=dmboundsuni(de,x,l,u);
      % left boundary of interval
    if (fl>0 and not((b1:=car b)<0))
     or(fl<0 and not((b1:=cdr b)>0))
     or not((pt:=l-fl/b1)<u)  then goto loop;
    if pt=l then goto precerr;
    q:=evaluateuni(e,x,pt);
    if not all and q*fl<0 then return regula!-falsi(e,x,l,fl,pt,q,acc);
    if abs q<acc then goto root_found;
    l:=pt; fl:=q;
      % right boundary of interval
    if (fu>0 and not((b1:=cdr b)>0))
     or(fu<0 and not((b1:=car b)<0))
     or not((pt:=u-fu/b1)>l) then goto loop;
    if pt=u then goto precerr;
    q:=evaluateuni(e,x,pt);
    if not all and q*fu<0 then return regula!-falsi(e,x,pt,q,u,fu,acc);
    if abs q<acc then goto root_found;
    u:=pt; fu:=q;
      % new point    
    pt :=(l+u)/2;  %midpoint
    q:=evaluateuni(e,x,pt);
    if not all and q*fu<0 then return regula!-falsi(e,x,l,fl,pt,q,acc);
       % generate new intervals
    if not(abs q<acc) then
     <<il2 :=find!-root!-addinterval(pt.q,u.fu,
               find!-root!-addinterval(l.fl,pt.q,il2));
       goto loop;
     >>;
  root_found:
    r:=pt.r; if not all then goto done;
    rr:=find!-root!-range(e,x,pt,acc);
    il2 :=find!-root!-addinterval(cdr rr,u.fu ,
               find!-root!-addinterval(l.fl,car rr,il2));
    goto loop;
  bottom:
    il1 :=il2;
    goto main_loop;
  done:
    r:=sort(r,function lessp);
    return make!-rdsolve!-result1(x,r);
  precerr:
    rederr "precision not sufficient for finding all roots";
  end;

symbolic procedure rd!-solve!-hardcase!-print(lev,il1);
 if !*trnumeric then
 << prin2t {"recursion level:",lev,"remaining intervals:",length il1};
    writepri( mkquote( 'list.
           for each i in il1 collect
             {'list,{'!*interval!*,prepf caar i,prepf cadr i},
                               dm!: prepf(cadr i-caar i)}),
            'only);
 >>;

symbolic procedure find!-root!-addinterval(l,h,il);
 dm!:  <<if car l < car h then il:=(l.h).il; il>>;

symbolic procedure find!-root!-range(e,x,p,acc);
  % p is a point in (l .. u) where abs e(x)<acc.
  % Find the next interval where e(x) > acc.
dm!:<<
     while not(abs(fl:=evaluateuni(e,x,pl:=pl-acc/2))>acc) do;
     while not(abs(fu:=evaluateuni(e,x,pu:=pu+acc/2))>acc) do;
     (pl.fl).(pu.fu)
    >> where pl=p,pu=p,fl=nil,fu:=nil;

if errorp errorset('(!:difference nil nil),nil,nil) then
<<
symbolic procedure !:difference(u,v);
   if null u then !:minus v else if null v then u
    else if u=v then nil
    else if atom u and atom v then u-v
    else dcombine(u,v,'difference);

symbolic procedure !:quotient(u,v);
   if null u or u=0 then nil
    else if null v or v=0 then rerror(poly,12,"Zero divisor")
    else if atom u and atom v
     % We might also check that remainder is zero in integer case.
     then if null dmode!* then u/v else dcombine(u,!:recip v,'times)
    else dcombine(u,v,'quotient);

>>;

endmodule;


module newton;  % root finding with generalized Newton methods.

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1992, all rights resrved
% Version 2: termination and damping criterion modified to affine
%            invariance.

fluid '(!*noequiv accuracy!*);
global '(iterations!* !*trnumeric erfg!*);


symbolic procedure rdnewton0(e,vars,p,n);
    %preparations for Newton iteration.
 (begin scalar l,jac,x,y,oldmode,q,!*noequiv;
    integer prec;
    if not memq(dmode!*,'(!:rd!: !:cr!:))then
       <<oldmode:=t; setdmode('rounded,t)>>;
    prec:=precision 0;
    p := for each x in p collect
         force!-to!-dm numr simp x;
    if !*trnumeric then lprim "computing symbolic Jacobian";
    eval list('matrix,mkquote list list('jacobian,n,n));
    for i:=1:n do for j:=1:n do
        setmatelem(list('jacobian,i,j),
         reval list('df,nth(e,i),nth(vars,j)));
    if !*trnumeric then lprim "inverting symbolic Jacobian";
    jac:= cdr reval'(quotient 1 jacobian);
    jac := for each r in jac collect
       for each c in r collect reval c;
    !*noequiv:=t;
    x := rdnewton1(e,jac,vars,p,'root);
    if oldmode then setdmode('rounded,nil);
    precision prec;
    if null x then rederr "no solution found";
    return 'list.
       for each p in pair(vars,x) collect
          list('equal,car p,cdr p);
  end) where !*roundbf=!*roundbf;

symbolic procedure rdnewton1(f,jac,vars,x,mode);
  begin scalar r,acc;
        integer n;
     if !*trnumeric then lprim "starting Newton iteration";
     acc := !:!:quotient(1,expt(10,accuracy!*));
     r:= rdnewton2(f,jac,vars,acc,x,mode,nil,nil);
     r:= for each x in r collect prepf x;
     return r;
 end;

symbolic procedure rdnewton2(f,jac,vars,acc,x,mode,low,high);
  % Algorithm for finding the root function system f
  % with technique of adaptively damped Newton.
  % f:     function to minimize (list of algebraic exprs);
  % jac:   symbolic inverted Jacobian;                    
  % vars:  variables (list of id's);
  % acc:   requested accuracy (e.g. 0.0000001)
  % x:     starting point (list of domain elements).
 dm!:
  begin scalar n,n0,n1,n2,e0,e1,dx,dx2,x1,g,dmp,delta,h;
        scalar dxold,dx2old,dmpold;
    integer count;
    e0 := list!-evaluate(f,vars,x);
  loop:
    count:=add1 count;
    if count>iterations!* then
    <<lprim "requested accuracy not reached within iteration limit";
      goto ready>>;

      % evaluate inverse Jacobian.
    g := matrix!-evaluate(jac,vars,x);
    dx := mat!*list(g,e0);    % the newton step.
    n0 := normlist dx;
    dmp:= 1;

  step:
      % evaluate function at new point.
    x1 := list!-list(x, scal!*list(dmp,dx));
    e1 := errorset({'list!-evaluate,mkquote f,
              mkquote vars,mkquote x1},nil,nil);
    if errorp e1 then goto contract else e1 := car e1;

        % anticipated next (simplified) Newton step.
    n1 := normlist(dx2:=mat!*list(g,e1));
    if n1=0 or n1<n0 then goto accept;
    if null dmpold then goto contract;

       % predict optimal damping factor
    h:= dmpold *(normlist list!-list(dx2old,dx) * normlist(dx))
        / (normlist(dxold) * normlist(dx2old));
    if h>1 then
    << dmp:=if h<10 then 1/h else 1/10;
       dmpold := nil; goto step>>;

  contract:
    dmp:=dmp/2;
    if dmp<acc then rederr "Newton method does not converge";
    goto step;

  accept:
    delta:=dmp*n0;
    x := x1; e0:=e1; n0:=n1;
    if low and high and (low>car x or high<car x) then return nil;
    dmpold := dmp; dxold := dx; dx2old := dx2;
    rdnewtonprintpoint(count,x,delta,e0);
    if n1>acc or dmp<1 then<<update!-precision(delta.e0); goto loop>>;

  ready:
    x := list!-list(x,dx2);
    return x;
  end;

symbolic procedure rdnewtonprintpoint(count,x,dx,e0);
  if !*trnumeric then
   begin integer n;
    writepri(count,'first);
    writepri(". residue=",nil);
    printsflist e0;
    writepri(", step length=",nil);
    writepri(mkquote prepf dx,'last);
    writepri(" at ",nil);
    printsflist x;
    writepri(" ",'last);
   end;

endmodule;


module steepstd;  % Optimization and root finding with method of
                  % steepest descent.


% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1991, all rights resrved

fluid '(!*noequiv accuracy!*);
global '(iterations!* !*trnumeric);

symbolic procedure rdmineval u;
  begin scalar e,l,vars,x,y,z,oldmode,p,!*noequiv;
    oldmode:=steepdecedmode(nil,nil);
    u := for each x in u collect reval x;
    u := accuracycontrol(u,6,40);
    e :=car u; u :=cdr u;
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if eqcar(x,'equal) then
       <<z:=cadr x; y:=caddr x>> else <<z:=x; y:=random 100>>;
       vars:=z . vars; p := y . p>>;
    vars := reversip vars; p := reversip p;
    x := steepdeceval1(e,vars,p,'min);
    steepdecedmode(t,oldmode);
    return list('list, car x, 'list.
       for each p in pair(vars,cdr x) collect
          list('equal,car p,cdr p));
  end;

put('num_min,'psopfn,'rdmineval);

symbolic procedure rdsolvestdeval u;
   % solving a system of equations via minimizing the
   % sum of sqares.
  begin scalar e,l,vars,x,y,z,oldmode,p,q,!*noequiv;
    oldmode:=steepdecedmode(nil,nil);
    u := for each x in u collect reval x;
    e :=car u; u :=cdr u;
      % construct the function to minimize.
    e:=if eqcar(e,'list) then cdr e else list e;
    q := 'plus . for each f in e collect
       list('expt,if eqexpr f then !*eqn2a f else f,2);
    q := prepsq simp q;
      % starting point & variables.
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if eqcar(x,'equal) then
       <<z:=cadr x; y:=caddr x>> else <<z:=x; y:=random 100>>;
       vars:=z . vars; p := y . p>>;
    vars := reversip vars; p := reversip p;
    x := steepdeceval1(q,vars,p,'root);
    steepdecedmode(t,oldmode);
    if null x then rederr "no solution found";
    return 'list.
       for each p in pair(vars,cdr x) collect
          list('equal,car p,cdr p);
  end;

symbolic procedure steepdecedmode(m,oldmode);
  if not m then % initial call
   << if !*complex then rederr
      "steepest descent method not applicable under complex";
    if not (dmode!*='!:rd!:)then
       <<oldmode:=t; setdmode('rounded,t)>>;
    {oldmode,precision 0,!*roundbf}>>
  else % reset to old dmode.
    <<if car oldmode then setdmode('rounded,nil);
      precision cadr oldmode;
      !*roundbf:=caddr oldmode>>;

symbolic procedure steepdeceval1(f,vars,x,mode);
  begin scalar e,g,r,acc;
        integer n;
     n := length vars;
    % establish the target function and the gradient.
     e := prepsq simp f;
     if !*trnumeric then lprim "computing symbolic gradient";
     g := for each v in vars collect prepsq simp list('df,f,v);
     !*noequiv:=t;
     if !*trnumeric then
      lprim "starting Fletcher Reeves steepest descent iteration";
     acc := !:!:quotient(1,expt(10,accuracy!*));
     x := for each u in x collect force!-to!-dm numr simp u;
     r:= steepdec2(e,g,vars,acc,x,mode);
     r:= for each x in r collect prepf x;
     return r;
 end;

symbolic procedure steepdec2(f,grad,vars,acc,x,mode);
  % Algorithm for finding the minimum of function f
  % with technique of steepest descent, version Fletcher-Reeves.
  % f:     function to minimize (standard quotient);
  % grad:  symbolic gradient (list of standard quotients);
  % vars:  variables (list of id's);
  % acc:   target precision (e.g. 0.000001)
  % x:     starting point (list of numerical standard forms).
  % mode:  minimum / roots type of operation
 dm!:
  begin scalar e0,e00,e1,e2,a,a1,a1a1,a2,a2a2,x1,x2,g,h,dx,w,
      delta,limit,gold,gnorm,goldnorm,multi;
    integer count,k,n;
    n:=length grad; multi:=n>1; n:=10*n;
    e00 := e0 := evaluate(f,vars,x);
    a1 := 1;
  init:
    k:=0;
      % evaluate gradient (negative).
    g := for each v in grad collect -evaluate(v,vars,x);
    gnorm := normlist g; h:=g;
  loop:
    count := add1 count; k:=add1 k;
    if count>iterations!* then
    <<lprim "requested accuracy not reached within iteration limit";
      % mode := nil;
      goto ready>>;

      % quadratic interpolation in direction of h in order
      % to find the minimum value of f in that direction.
    a2 := nil;
   l1:
    x1 := list!+list(x, scal!*list(a1,h));
    e1 := evaluate(f,vars,x1);
    if e1>e0 then <<a2:=a1; x2:=x1; e2:=e1;
            a1 := a1/2; goto l1>>;
    if a2 then goto alph;
   l2:
    a2 := a1+a1;
    x2 := list!+list(x, scal!*list(a2,h));
    e2 := evaluate(f,vars,x2);
    if e2<e1 then <<a1:=a2; e1:=e2; goto l2>>;
   alph:   %compute lowest point of parabel
    if abs(e1-e2)<acc then goto ready; % constant?
    a1a1:=a1*a1; a2a2:=a2*a2;
    a:= (a1a1-a2a2)*e0 + a2a2*e1 - a1a1*e2 ;
    a:= a/((a1-a2)*e0 + a2*e1-a1*e2);
    a:= a/2;

     % new point
    dx:=scal!*list(a,h); x:=list!+list(x, dx);
    e0 := evaluate(f,vars,x);
    if e0>e1 then % a1 was better than interpolation.
    <<dx:=scal!*list(a1-a ,h); x:=list!+list(x,dx);
      e0 := e1; dx:=scal!*list(a1,h)>>;
    delta:= normlist dx;
    steepdecprintpoint(count,x,delta,e0,gnorm);
    update!-precision list(delta,e0,gnorm);

      % compute next direction;
    if k>n then goto init; % reinitialize time to time.
    gold := g; goldnorm:= gnorm;
    g := for each v in grad collect -evaluate(v,vars,x);
    gnorm := normlist g;
    if gnorm<limit then goto ready;
    h := if not multi then g else
        list!+list(g,scal!*list(gnorm/goldnorm,h));
    if mode='min and gnorm > acc or
       mode='root and e0 > acc then goto loop;

  ready:
    if mode='root and not(abs e0 < acc) then
     <<lprim "probably fallen into local minimum of f^2";
        return nil>>;
    return e0 . x;
  end;

symbolic procedure steepdecprintpoint(count,x,d,e0,ng);
 if !*trnumeric then
   begin integer n;
    writepri(count,'first);
    writepri(". residue=",nil);
    writepri(mkquote prepf e0,nil);
    writepri(", gradient length=",nil);
    writepri(mkquote prepf ng,'last);
    writepri(" at (",nil);
    for each y in x do
       <<if n>0 then writepri(" , ",nil);
         n:=n+1;
         writepri(mkquote prepf y,nil)>>;
    writepri(")",nil);
    writepri(", stepsize=",nil);
    writepri(mkquote prepf d,nil);
    writepri("",'last);
   end;

endmodule;


module bounds;   % Upper and lower bound of a 
% multivariate functional expression in a n-dimensional interval.

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1991, all rights resrved

put('bounds,'psopfn,'boundseval);
put('bounds!-rd,'psopfn,'boundsevalrd);
put('bounds,'numericfn,'bounds!-rd);

fluid '(boundsdb!* !*msg boundspolynomialdb!*
       !*boundsfullcheck boundsknown!*);

symbolic procedure boundsevalrd u;
   begin scalar dmode!*;
     setdmode('rounded,t);
     return boundseval u;
   end;

symbolic procedure boundseval u;
 (begin scalar db,y,l,f;
    u := for each x in u collect reval x;
    f := car u; u := cdr u;
    if u and eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
    <<if not eqcar(x,'equal) then typerr(x,"domain description");
      y := revalnuminterval(caddr x,nil);
      l:=list(simp car y,simp cadr y);
      db := (cadr x . minsq l . maxsq l).db>>;
    u := boundseval1(f,db);
    return mkinterval(prepsq car u,prepsq cdr u);
  end) where !*roundbf=!*roundbf;

symbolic procedure boundserr(a,b);
    if not !*msg then error(99,'bounds) else
    if a then typerr(a,b) else rederr b;

symbolic procedure boundseval0 (f,db,!*msg);
  % internal calling convention:
  %  f            algebraic expression,
  %  db           list of (x . lo . hi), where
  %               lo and hi represent the current interval for
  %               variable x as standard quotients.
  %  result is a pair of standard quotients representing
  %               the bounds for f.
    boundseval1(f,db);

symbolic procedure boundseval1(u,boundsdb!*);
    <<if u member boundsknown!* then
      !*boundsfullcheck := nil else
      !*boundsfullcheck := boundsknown!* := u.boundsknown!*;
      boundseval2 u>>;

symbolic procedure boundseval2 u;
  % Main driver for bounds evaluation
  if adomainp u then <<u:=simp u;u.u>> else
  boundspolynomial u or
  begin scalar v,w,fcn;
    if (v:=assoc(u,boundsdb!*))then return cdr v;
    if idp u and (fcn:=get(u,dmode!*)) then
      <<v:=apply(fcn,nil); v:=(v ./1).(v ./1); goto done>>;
    if atom u then goto err;
    if car u='expt and fixp caddr u and caddr u>0 then
    <<v:=boundsexpt!-int(boundseval2 cadr u,caddr u);
      goto done>>;
    w := for each y in cdr u collect boundseval2 y;
    v:= if cdr w and cddr w and (fcn:=get(car u,'boundsevaln))
           then apply1(fcn,w)
         else if length w>2 then t
         else if (fcn:=get(car u,'boundseval1))
          then apply2(fcn,car u,car w)
         else if (fcn:=get(car u,'boundseval))
          then if null cdr w then apply1(fcn,car w)
                else apply2(fcn,car w,cadr w)
         else if cdr w then t
         else boundsoperator(car u,car w);
    if v=t then goto err2;
    if null v or
       not domainp numr car v or not domainp denr car v or
       not domainp numr cdr v or not domainp denr cdr v then goto err;
     % cache result for later usage.
  done:
    boundsdb!*:= (u.v).boundsdb!*;
    return v;
err: boundserr(nil,"unbounded in range");
err2: boundserr(nil,"bounds confusion");
  end;

symbolic procedure boundsoperator(fcn,u);
  % general external interface: function flagged decreasing,
  % increasing of explicit bounds given.
   if flagp(fcn,'increasing) then boundsincreasingfn(fcn,u) else
   if flagp(fcn,'decreasing) then boundsdecreasingfn(fcn,u) else
   if get(fcn,'upperbound) and get(fcn,'lowerbound) then
       simp get(fcn,'lowerbound) . simp get(fcn,'upperbound)
    else nil;

symbolic procedure boundsplus u;
   if null cdr u then car u else
   boundsplus2(car u,boundsplus cdr u);

symbolic procedure boundsplus2(u,v);
    addsq(car u,car v) . addsq(cdr u,cdr v);

put('plus,'boundsevaln,'boundsplus);
put('plus,'boundseval,'boundsplus2);

symbolic procedure boundsdifference(x,y);
   subtrsq(car x,cdr y).subtrsq(cdr x,car y);

put('difference,'boundseval,'boundsdifference);

symbolic procedure boundsminus(u);
   negsq cdr u . negsq car u;

put('minus,'boundseval,'boundsminus);

symbolic procedure boundstimes u;
   if null cdr u then car u else
    boundstimes2(car u,boundstimes cdr u);

symbolic procedure boundstimes2(x,y);
  begin scalar z;
      % first handle simple cases
   if not minusf numr car x and not minusf numr car y then
        return multsq(car x,car y) . multsq(cdr x,cdr y);
   if minusf numr cdr x and minusf numr cdr y then
        return multsq(cdr x,cdr y) . multsq(car x,car y);
   z:=list(multsq(car x,car y), multsq(car x,cdr y),
           multsq(cdr x,car y), multsq(cdr x,cdr y));
   return minsq z . maxsq z;
  end;

symbolic procedure minsq l;
  begin scalar x,z;
   x := car l;
   for each y in cdr l do
    if minusf numr subtrsq(y,x) then x:=y;
   return x;
  end;

symbolic procedure maxsq l;
  begin scalar x,z;
   x:= car l;
   for each y in cdr l do
    if minusf numr subtrsq(x,y) then x:=y;
   return x;
  end;

symbolic procedure sqgreaterp(x,y);
     minusf numr subtrsq(y,x);

put('times,'boundsevaln,'boundstimes);
put('times,'boundseval,'boundstimes2);

symbolic procedure boundsexp u;
   boundsexpt(b . b,u) where b = simp 'e;

put('exp,'boundseval,'boundsexp);

symbolic procedure boundsexpt!-int(b,ex);
  % Form x ** n. Look for even exponents.
  begin scalar hi,lo,bh,bl;
    bl := exptsq(car b,ex); bh := exptsq(cdr b,ex);
    if not evenp ex then return bl.bh;
    lo := minusf numr cdr b;
    hi := not minusf numr car b;
    return
        if hi then bl.bh else  % both had been positive
        if lo then bh.bl else  % both had been negative: invert
             (nil ./ 1) . maxsq list(bl,bh);
   end;

symbolic procedure boundsexpt!-const(b,ex);
  % Form x ** constant, including fractional exponents.
  begin scalar l;
    l := list(bounds!-evalfcn2('expt,car b,ex),
              bounds!-evalfcn2('expt,cdr b,ex));
    return minsq l . maxsq l;
  end;

symbolic procedure boundsexpt(b,e);
   if car e=cdr e then boundsexpt!-const(b,car e) else
  % Form x ** y. x>0 only.
  % Either monotonous growing or falling: pick extremal points.
  begin scalar l;
    if sqgreaterp(nil./1,car b) then return nil;
    l := list(bounds!-evalfcn2('expt,car b,car e),
              bounds!-evalfcn2('expt,car b,cdr e),
              bounds!-evalfcn2('expt,cdr b,car e),
              bounds!-evalfcn2('expt,cdr b,cdr e));
    return minsq l . maxsq l;
  end;

symbolic procedure boundsprepsq u; prepsq car u . prepsq cdr u;

put('expt,'boundseval,'boundsexpt);

symbolic procedure boundsquotient(n,d);
  begin scalar l;
    if minusf numr multsq(car d,cdr d) then
      boundserr(nil,"unbounded in range");
    l := list(quotsq(car n,car d),quotsq(car n,cdr d),
              quotsq(cdr n,car d),quotsq(cdr n,cdr d));
    return minsq l . maxsq l;
  end;

put('quotient,'boundseval,'boundsquotient);

symbolic procedure boundsabs u;
 if evalgreaterp(prepsq car u,0) then u else
 if evalgreaterp(0,prepsq cdr u) then negsq cdr u . negsq car u else
  (nil ./1) . maxsq list (negsq car u,cdr u);

put('abs,'boundseval,'boundsabs);

symbolic procedure boundsincreasingfn(fn,u);
  % Bounds for an increasing function. Out-of -range problems will
  % be detected either by the function evaluation or finally by
  % the main driver if the result is not an interval with numeric
  % bounds.
  bounds!-evalfcn1(fn,car u) . bounds!-evalfcn1(fn,cdr u);

put('log,'boundseval1,'boundsincreasingfn);
put('asin,'boundseval1,'boundsincreasingfn);
put('atan,'boundseval1,'boundsincreasingfn);
put('sqrt,'boundseval1,'boundsincreasingfn);

symbolic procedure boundsdecreasingfn(fn,u);
       boundsincreasingfn(fn,cdr u.car u);

put('acos,'boundseval1,'boundsdecreasingfn);
put('acot,'boundseval1,'boundsdecreasingfn);

symbolic procedure boundssincos(fcn,n);
   % Reason if one of the turn points (n*pi) is in the
   % range. If yes, include the corresponding value 1 or -1.
   % Otherwise compute the range spanned by the end points.
  begin scalar lo,hi,alo,ahi,!1pi,!2pi,!3pi,l,m;
        integer s;
     lo := car n; hi := cdr n;
     <<setdmode('rounded,t);
       !1pi := simp 'pi;
       setdmode('rounded,nil);
     >> where dmode!* =nil;
     if not domainp numr !1pi then goto trivial;
     !2pi := addsq(!1pi,!1pi); !3pi := addsq(!1pi,!2pi);
       % convert sin to cos
     if fcn = 'sin then <<m :=multsq(!1pi, -1 ./ 2);
         lo := addsq(lo,m); hi := addsq(hi,m)>>;
     m := negsq multsq(!2pi,fixsq quotsq(lo,!2pi));
       % move interval to main interval
     lo:=addsq(lo,m); hi:=addsq(hi,m);
     if minusf numr lo then
     <<lo := addsq(lo,!2pi); hi := addsq(hi,!2pi)>>;
     if null numr lo or sqgreaterp(hi,!2pi) then l:= (1 ./ 1) . l;
     if (sqgreaterp(!1pi,lo) and  sqgreaterp(hi,!1pi))
      or(sqgreaterp(!3pi,lo) and  sqgreaterp(hi,!3pi))
       then l := (-1 ./ 1) . l;
     if l and cdr l then goto trivial;
     l:=bounds!-evalfcn1('cos,lo) .bounds!-evalfcn1('cos,hi).l;
     if smemq('cos,l) then goto trivial;
     return minsq l . maxsq l;
 trivial:
     return ((-1)./1) . (1 ./ 1);
  end;

symbolic procedure fixsq u;
         bounds!-evalfcn1('fix,u);

symbolic procedure bounds!-evalfcn1(fcn,u);
   (if domainp numr u and denr u=1
         and (x:=valuechk(fcn,list numr u))
     then x else simp list(fcn,prepsq u)) where x=nil;

symbolic procedure bounds!-evalfcn2(fcn,u,v);
   (if domainp numr u and denr u=1 and
       domainp numr v and denr v=1
         and (x:=valuechk(fcn,list(numr u,numr v)))
     then x else simp list(fcn,prepsq u)) where x=nil;

symbolic procedure agreaterp(u,v);
   (lambda x;
    if not atom denr x or not domainp numr x
      then error(99,"number")
     else numr x and !:minusp numr x)
        simp!* list('difference,v,u);

symbolic procedure boundssin u;
       boundssincos('sin,u);

symbolic procedure boundscos u;
       boundssincos('cos,u);

put('sin,'boundseval,'boundssin);
put('cos,'boundseval,'boundscos);

symbolic procedure boundstancot(fcn,n);
  begin scalar lo,hi,alo,ahi,x,m,ppi; integer s;
     lo := car n; hi := cdr n;
     ppi := rdpi!*() ./ 1;
     if sqgreaterp(subtrsq(lo,hi),ppi) then goto no;
     lo:=bounds!-evalfcn1(fcn,lo);
     hi:=bounds!-evalfcn1(fcn,hi);
     if fcn='cot then <<x:=lo;lo:=hi;hi:=x>>;
     if not sqgreaterp(lo,hi) then return lo.hi;
no:  return boundserr(nil,"unbounded in range");
  end;


symbolic procedure boundstan u;
    boundstancot('tan,u);

symbolic procedure boundscot u;
    boundstancot('cot,u);

put('tan,'boundseval,'boundstan);
put('cot,'boundseval,'boundscot);

symbolic procedure boundsmaxmin u;
  if null cdr u then car u else
       boundsmaxmin2(car u,boundsmaxmin cdr u);

symbolic procedure boundsmaxmin2(x,y);
  (if sqgreaterp(car x,car y) then car y else car x).
  (if sqgreaterp(cdr x,cdr y) then cdr x else cdr y);

put('max,'boundsevaln,'boundsmaxmin);
put('min,'boundsevaln,'boundsmaxmin);
put('max,'boundseval,'boundsmaxmin2);
put('min,'boundseval,'boundsmaxmin2);

% for univariate polynomials we look at the extremal points and the
% interval ends. Assuming that the same expression will be investigated
% with different intervals, we store the knowledge about the polynomial.

symbolic procedure boundspolynomial e;
 dm!: begin scalar x,u,z1,z2,lo,hi,mx,mn,d,v,l;
    if dmode!* neq '!:rd!: then return nil;
    if(u:=assoc(e,boundspolynomialdb!*)) then goto old;
    if null !*boundsfullcheck or
       null(x:=boundspolynomialtest e) then return nil;
    d := realroots{reval{'df,e,x}};
    u := x. for each s in cdr d collect
      <<d:=numr simp caddr s;
        d . evaluate(e,{x},{d})>>;
    u:=e.u;
    boundspolynomialdb!*:=u.boundspolynomialdb!*;
 old:
    x:=cadr u; u :=cddr u;
    v := assoc(x,boundsdb!*);
    if null v then return nil;
    lo:=numr cadr v; hi:=numr cddr v;
    l:={evaluate(e,{x},{lo})./1 , evaluate(e,{x},{hi}) ./1};
    for each p in u do
      if lo<car p and car p<hi then l:= (cdr p./1).l;
    return minsq l . maxsq l;
  end;

symbolic procedure boundspolynomialtest e;
 (pairp r and car r) where r=boundspolynomialtest1(e,nil);

symbolic procedure boundspolynomialtest1(e,x);
   if adomainp e then x or t else
   if atom e then boundspolynomialtestvariable(e,x) else
   if car e eq 'expt then
       fixp caddr e and
          boundspolynomialtestvariable(cadr e,x)
      else
   if car e eq 'minus then boundspolynomialtest1(cadr e,x)
      else
   if car e  eq 'plus or car e eq 'times then
    begin scalar r;
      e:=cdr e; r:=t;
      while e and r do
      <<r:=boundspolynomialtest1(car e,x);
        if r neq t then x:=r; e:=cdr e>>;
      return x or r;
    end
      else boundspolynomialtestvariable(e,x);

symbolic procedure boundspolynomialtestvariable(e,x);
  if x and e=car x then x else
  if x then nil else
  if assoc(e,boundsdb!*) then {e};

%--------------------------------------------------------------
%
% support for compilation
%
%--------------------------------------------------------------

fluid '(boundsdb!* boundscompv!* boundscompcode!*);

symbolic procedure boundscompile(e,v);
  % compile bounds evaluation function for expression e,
  % input intervals v.
  begin scalar boundsdb!*,boundscompv!*,boundscompcode!*,u;
    boundsdb!*:=for each x in v collect x.x;
    u:=boundscompile1(e,nil);
    return
     {'lambda,v,
      'prog . boundscompv!* .
         append(reversip boundscompcode!*,{{'return,u}})};
  end;

symbolic procedure boundscompile1(u,flag);
  % Main driver for bounds compilation
  if adomainp u then <<u:=simp u;mkquote(u.u)>> else
  begin scalar v,w,fcn,var,q;
    if (v:=assoc(u,boundsdb!*))then return cdr v;
    if idp u and (fcn:=get(u,dmode!*)) then
      <<v:=apply(fcn,nil); v:=mkquote((v ./1).(v ./1)); goto done>>;
    if atom u then goto err1;
    if car u='expt and fixp caddr u and caddr u>0 then
    <<v:={'boundsexpt!-int,boundscompile1(cadr u,t),caddr u};
      goto done>>;
    w := for each y in cdr u collect boundscompile1(y,t);
    v:= if length w>2 and (fcn:=get(car u,'boundsevaln))
              then {fcn,'list.w} else
        if length w>2 then t else
        if (fcn:=get(car u,'boundseval1))
               then {fcn,mkquote car u,car w} else
        if (fcn:=get(car u,'boundseval))
               then if null cdr w then {fcn,car w}
                     else {fcn,car w,cadr w}
         else if cdr w then t
         else {'boundsoperator,car u,car w};
  done:
    if v=t then goto err2;
    if not flag then return v;
    var:=gensym(); boundscompv!* := var.boundscompv!*;
    boundscompcode!*:={'setq,var,v}.boundscompcode!*;
     % cache result for later usage.
    boundsdb!*:= (u.var).boundsdb!*;
    return var;
err1: typerr(u,"bounded value");
err2: typerr(u,"expression to be compiled for bounds");
  end;

endmodule;


module numint; % compute a numerical integral

% Author: H. Melenk, 

     % Version March 1993

     % strategy for univariate integrals:
     %
     % 1. look for exact antiderivative, use that if it is 
     %          bounded in the interval,
     %
     % 2. try Chebyshev fit,
     %
     % 3. use M. Wulkow's adaptive multilevel quadrature.

fluid '(!*noequiv accuracy!* singularities!*);
global '(iterations!* !*trnumeric);

symbolic procedure intrdeval u;
     % interface function;
  begin scalar e,l,vars,y,p,singularities!*;
        integer n;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,3,20);
    e :=car u; u :=cdr u;
      % variables and interval's.
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if not eqcar(x,'equal) then typerr(x,"interval bounds");
       y:=revalnuminterval(caddr x,t);
       vars:=cadr x.vars;
       p:=(cadr x. car y. cadr y).p>>;
    return intrd1(e,vars,p);
  end;

put('num_int,'psopfn,'intrdeval);

symbolic procedure intrd1 (e,vars,p);
  if cdr vars then intrd2(e,vars,p) else
     % Univariate case.
   begin scalar fcn,x,u,w,acc,lo,hi,uu,loo,hii,oldmode,cbound;
     integer ord;
    fcn := e;
    x := car vars;
    lo := car cdar p; hi := cdr cdar p;
    if evalgreaterp(hi,lo) then <<loo:=lo; hii:=hi>>
      else <<loo:=hi; hii:=loo>>;
     % first attempt: look for bounded antiderivative.
    u := reval{'int,reval fcn,x} where dmode!*=nil,!*rounded=nil;
    oldmode:=switch!-mode!-rd nil;
    w:= errorset(
        {'boundseval,mkquote{u, {'equal,x,{'!*interval!*, loo, hii}}}},
                             nil,nil) where !*msg=nil;
    switch!-mode!-rd oldmode;

    if not smemq('int,u) and not errorp w then
      <<
        if !*trnumeric then prin2t "Using bounded antiderivative";
        oldmode:=switch!-mode!-rd nil;
        u:=simp u;
        w:= prepsq subtrsq(subsq(u,{x . hi}),subsq(u,{x . lo}));
        switch!-mode!-rd oldmode;
        return w;
      >>;
    if !*trnumeric then prin2t "No bounded antiderivative found";

    % second attempt: look for a Chebyshev fit.
    w := nil;
    oldmode:=switch!-mode!-rd nil;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    cbound := dm!: (acc /(hii - loo));
    ord := 20;
  chebloop:
    u := chebcoeffs(fcn,x,loo,hii,ord);
    uu := reverse u;
    if int!-chebconverges(uu,acc,cbound) then
    <<
      u:=chebint(u,nil,loo,hii);
      w:=aeval{'difference,chebeval(u,nil,loo,hii,hi),
                           chebeval(u,nil,loo,hii,lo)};
    >>;
    if null w and ord<60 then <<ord:=ord+20; goto chebloop>>;
    switch!-mode!-rd oldmode;
    if w then
    <<if !*trnumeric then
       << prin2 "Using Chebyshev approximation or order ";
          prin2t ord>>;
      return w;
    >>;
    if !*trnumeric then
    <<
       prin2t "No usable Chebyshev approximation found";
       prin2t "Starting adaptive multilevel quadrature";
    >>;
    return intrd2 (e,vars,p);
  end;

symbolic procedure int!-chebconverges(u,acc,ab);
  % Test the (reverse) Chebyshev coefficients u for convergence. 
  % acc: current accuracy (epsilon), ab: absolute upper bound for
  % trailing coefficient(s).
  begin scalar mx;
    mx := int!-chebmax(u,nil);
    dm!:(mx := mx*acc);
    return dm!:(abs car u < mx
        and abs cadr u < mx and abs caddr u < mx
        and abs car u < ab and abs cadr u < ab
        and abs caddr u < ab);
  end;

symbolic procedure  int!-chebmax(u,mx);
   if null u then mx else
   int!-chebmax(cdr u,if dm!:(w>mx) then w else mx)
          where w=dm!: abs car u;


% ----------------- adaptive multilevel quadrature

symbolic procedure intrd2 (e,vars,p);
   begin scalar acc,r,oldmode,!*noequiv;
    oldmode:=switch!-mode!-rd nil;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    e := reval e;
    r := if null cdr p then intrduni(e,p,acc) else
          intrdmulti(e,p,acc);
    switch!-mode!-rd oldmode;
    return r;
  end;

symbolic procedure intevaluate1(e,x,v);
  (if atom a then <<singularities!*:=v.singularities!*; 0>>
     else car a)where a=evaluate!*(e,list x,list v);

symbolic procedure intevaluate(e,x,v);
  (if atom a then <<singularities!*:=v.singularities!*; 0>>
     else car a)where a=evaluate!*(e,x,v);

symbolic procedure intrduni (e,p,acc);
   % univariate integral.
 dm!:
  begin scalar lo,flo,hi,fhi,x,d,u,eps,i,il,r,int,oldint;
     integer n,k;
    x := car p; p:= cdr p;
   % evaluate the interval.
    lo := cadr x; hi := cddr x; x := car x;
    lo:=force!-to!-dm lo; hi:=force!-to!-dm hi;
    if hi=lo then return force!-to!-dm nil;
      % initial interval;
    il := list intrdmkinterval(e,x,
               (lo.intevaluate1(e,x,lo)),
               (hi.intevaluate1(e,x,hi)),1);
    int:=car lastpair car il;
  loop:
    k:=add1 k; n:=add1 n;
    if remainder(n,25)=0 then intrdmess(n,il);
   % divide the interval with the biggest local error;
    i:=car il; il:=cdr il; u:=intrdintersect(e,x,i);
    if null u then
        <<il:=i.il;
          intrdabortmsg(list car cadr i,list x,intcollecteps il);
           goto ready>>;
    for each q in u do il := intrdsortin(q,il);
    if k<5 then goto loop;

   % control variation of result
    if n=5 then
    <<int:=intcollectint il;acc:=acc*abs int>>; % adjust accuracy
    k:=0; eps := intcollecteps il;
    if eps>acc then goto loop;
  ready:
    return intcollectint il;
  end;

symbolic procedure intrdabortmsg(pt,vars,eps);
   <<writepri(
"requested accuracy cannot be reached within iteration limit",'only);
     writepri("critical area of function near to ",'first);
     writepri(mkquote('list . for each u in pair(vars,pt)
       collect list('equal,car u,prepf cdr u)),'last);
     writepri("current error estimate: ",'first);
     writepri(mkquote prepf eps,'last);
   >>;

symbolic procedure intcollectint il;
 dm!: <<for each i in cdr il do r:= car lastpair i + r; r>>
          where r=car lastpair(car il);

symbolic procedure intcollecteps il;
 dm!: <<for each i in cdr il do r:= car i + r; r>>
          where r=car(car il);

symbolic procedure intrdsortin(u,l);
   % sort interval u into l such that l is sorted with descending car.
 dm!:( if null l then list u else
       if car u < caar l then car l . intrdsortin(u,cdr l) else
       u . l);

symbolic procedure intrdintersect(e,x,i);
  begin scalar d,plo,pmi,phi,lev;
    i:= cdr i;
    plo := car i; i := cdr i;
    pmi := car i; i := cdr i;
    phi := car i; i := cdr i;
    d   := car i; dm!:(d:=d/2);
    lev := cadr i +1;
    if lev>iterations!* then return nil;
    return list (intrdmkinterval(e,x,plo,pmi,lev) ,
                 intrdmkinterval(e,x,pmi,phi,lev) );
   end;

symbolic procedure intrdmkinterval(e,x,lo,hi,lev);
  dm!: begin scalar eps,it,is,mid,mi,flo,fhi,fmid,d,u;
    d := car hi- car lo;
    mid := (car hi+ car lo)/2;
    flo := cdr lo; fhi := cdr hi;
    fmid:=intevaluate1(e,x,mid);
    mi:=mid.fmid;
    if u:=intrdunisingular(lo,mi,hi) then
         <<flo:=car u; fmid:=cadr u; fhi:=caddr u>>;
    it := d*(flo+ 2fmid + fhi) / 4;         % order 1 integral;
    is := d*(4*flo + 16*fmid + 4*fhi)/24;   %simpson integral;
    eps:= abs(is-it);
                                            % interval record:
    return list (eps,                       % local error contrib
                 lo, mi, hi,                % the 3 points
                 d, lev, is);               % length and simpson.
  end;

symbolic procedure intrdunisingular(lo,mi,hi);
  % local extraploation for singular points.
  if singularities!* then
   begin scalar slo,smi,shi,fl,fm,fh,u;
    slo:=member(car lo,singularities!*);
    smi:=member(car mi,singularities!*);
    shi:=member(car hi,singularities!*);
    if not slo and not smi and not shi then return nil;
    if slo and smi and shi then rederr "too many singularities";
    fl:=cdr lo; fm:=cdr mi; fh:=cdr hi;
    dm!:( u := 2*(fl+fm+fh) );
    return list(if slo then u else fl,
                if smi then u else fm,
                if shi then u else fh);
  end;

symbolic procedure intrdmulti(e,p,acc);
  % multivariate integral.
 dm!:
  begin scalar vars,x,y,u,eps,i,il,r,int,oldint;
     integer n,k,dim;
    dim:=length p;
    il:=intrdmkinitcube(e,p);
    vars := car il; il:= cdr il;
  loop:
    k:=add1 k; n:=add1 n;
    if remainder(n,25)=0 then intrdmess(n,il);
   % divide the simplex with the biggest local error;
    i:=car il; il:=cdr il; u:=intrdrefinecube(e,vars,i);
    if null u then
        <<il:=i.il;
          intrdabortmsg(car cadr i,vars,intcollecteps il);
           goto ready>>;
    for each q in u do il := intrdsortin(q,il);
    if k<5 then goto loop;

   % control variation of result
    if n=5 then
    <<int:=intcollectint il;acc:=acc*abs int>>; % adjust accuracy
    k:=0; eps := intcollecteps il;
    if eps>acc then goto loop;
  ready:
    return intcollectint il;
  end;

symbolic procedure intrdmkinitcube(e,p);
  % make initial element.
   begin scalar vol,points,vars,il,lo,hi,x,y;
    vol:=1;
    for each z in p do
    <<vars:=car z.vars;
      x:=force!-to!-dm cadr z;
      y:=force!-to!-dm cddr z;
      lo:=x.lo; hi:=y.hi;
      vol:= dm!:( abs(x-y)*vol );
    >>;
    lo:=lo.intevaluate(e,vars,lo);
    hi:=hi.intevaluate(e,vars,hi);
    il:=list intrdmkcube(e,vars,lo,hi,vol,1);
    return vars.il;
  end;

symbolic procedure intrdrefinecube(e,vars,i);
 dm!:  % divide cube into 2.
  begin scalar plo,phi,lo,hi,nlo,nhi,vol,x,y,xnew;
     integer m,lev;
    i:=cdr i;
    plo:=car i; i:=cdr i; phi:=car i; i:=cdr i;
    vol:= car i / 2; lev := add1 cadr i;
    if lev > iterations!* then return nil;
    lo:=car plo; hi:=car phi;
      % select longest edge of interval;
    m := 1; x := car hi-car lo;
    for j:=2:length lo do
      if x<(y:=(nth(hi,j)-nth(lo,j))) then<<m:=j;x:=y>>;
    nlo := append(lo,nil); nhi := append(hi,nil);
    xnew:=(nth(hi,m) + nth(lo,m)) /2 ;
    nth(nlo,m):=xnew; nth(nhi,m):=xnew;
    nlo := nlo.intevaluate(e,vars,nlo);
    nhi := nhi.intevaluate(e,vars,nhi);
    return list(
      intrdmkcube(e,vars,plo,nhi,vol,lev),
      intrdmkcube(e,vars,nlo,phi,vol,lev));
  end;

symbolic procedure intrdmkcube(e,vars,lo,hi,vol,lev);
   % make a fresh cube record.
  dm!: begin scalar u,eps,it,is,mid,fmi,flo,fhi;
    flo:= cdr lo; fhi:= cdr hi;
    mid:=list!+list(car lo,car hi); mid:=scal!*list(1/2,mid);
    fmi:= intevaluate(e,vars,mid);
    if u:=intrdunisingular(lo,mid.fmi,hi) then
         <<flo:=car u; fmi:=cadr u; fhi:=caddr u>>;
    is:=flo+fhi;
    it:=is*vol/2; is:=(2*fmi+is)*vol/4;
    eps := abs(is-it);
    return list(eps,lo,hi,vol,lev,is);
  end;

symbolic procedure intrdmess(n,il);
  if !*trnumeric then
   <<writepri(2*n,'first); writepri(" intervals, integral=",nil);
     writepri(mkquote prepf intcollectint il,nil);
     writepri(", error estimate=",nil);
     writepri(mkquote prepf intcollecteps il,nil);
   >>;

symbolic procedure prinxt l;
  begin integer i;
    for each x in l do
     if not eqcar(x,'!:rd!:) then prin2 x else
     if atom cdr x then prin2 x else
       <<prin2 (float cadr x *expt(10.0,cddr x));
         tab (i:=20+i)>>;
     terpri();
   end;

endmodule;


module numfit; % approximation of functions with least spares.

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1991, all rights resrved

fluid '(!*noequiv accuracy!* singularities!*);
global '(iterations !*trnumeric);

symbolic procedure fiteval u;
     % interface function;
  begin scalar a,b,c,d,e,r,l,q,x,v,var,pars,fcn,fl,basis,pts,
        grad,oldmode,!*noequiv;
        integer n,i;
    if not(length u=3) then goto synerr;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,6,200);
    fcn :=car u; u :=cdr u;
    if eqcar(fcn,'list) then fl:=cdr fcn;
    basis:=car u; u :=cdr u;
    if not eqcar(basis,'list) then goto synerr;
    basis := for each x in cdr basis collect simp reval x;
    var:= car u; u := cdr u;
    if not eqcar(var,'equal) then goto synerr;
    if not eqcar(pts:=caddr var,'list) then goto synerr;
    var := cadr var; pts:=cdr pts; n:=length pts;
    if !*trnumeric then prin2t "build generic approximation function";
    a:=nil./1;
    for each b in basis do
    <<v:=gensym(); pars:=v.pars; a:=addsq(multsq(simp v,b),a)>>;
      % generate the error expression
    oldmode:=switch!-mode!-rd nil;!*noequiv:=nil;
    b:=a:=simp prepsq a;
    fcn:=simp if null fl then fcn else 'dummy; e:=nil./1;
    for each p in pts do
    <<i:=i+1;l:=list(var.p); % value to be substituted.
      if fl then l:=('dummy . reval nth(fl,i)).l;
      q:=subtrsq(subsq(fcn,l),subsq(b,l));
      e:=addsq(e,multsq(q,q))>>;
    e:=prepsq e;
    if !*trnumeric then
      <<lprim "error function is:";writepri(mkquote e,'only)>>;
      % find minimum.
      % build gradient.
    grad:='list . for each v in pars collect reval {'df,e,v};
    r:=rdsolveeval
       list (grad,'list . for each p in pars collect list('equal,p,0),
                  {'equal,'accuracy,accuracy!*},
                  {'equal,'iterations,iterations!*});
      % resubstitute into target function.
    if !*trnumeric then lprim "resubstituting in approximating form";
    l:=nil; pars:= nil;
    for each p in cdr r collect
       <<x:=caddr p; pars:=x.pars; l:=(cadr p.x).l>>;
    a:=prepsq subsq(a,l);
    switch!-mode!-rd oldmode;
    return list('list, a ,'list . pars);
  synerr:
    rederr "illegal parameters in fit";
  end;

put('num_fit,'psopfn,'fiteval);

endmodule;


module chebysh;  % Chebyshev approximation.
   % Literature: Press, Flannery, Teukolski, Vetterling: Numerical
   %                Recipes, Cambridge University Press.

% Author: Herbert Melenk

% Copyright (c) 1993 ZIB Berlin, RAND.  All rights reserved.

put('chebyshev_fit,'psopfn,'(lambda(u)(chebysheveval u 'fit)));
put('chebyshev_eval,'psopfn,'(lambda(u)(chebysheveval u 'eval)));
put('chebyshev_int,'psopfn,'(lambda(u)(chebysheveval u 'int)));
put('chebyshev_df,'psopfn,'(lambda(u)(chebysheveval u 'df)));

symbolic procedure chebysheveval(u,mode);
     % interface function;
  begin scalar w,x,c,e,l,y,p,r,oldmode;
        integer n;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,3,20);
    e :=car u; u :=cdr u;
      % variable and interval.
    w := car u;
    if not eqcar(w,'equal) then typerr(w,"interval bounds");
    oldmode:=switch!-mode!-rd nil;
    y:=revalnuminterval(caddr w,t);
    x:=cadr w;
    if mode = 'fit then
    <<
       n:=if cdr u then ieval cadr u else 20;
       c := chebcoeffs(e,x,car y,cadr y,n);
       r:= {'list,
            chebpol(c,x,car y,cadr y),
            'list . c}; %  for each q in c collect prepf q};
    >>
    else if mode = 'eval then
    <<
       if null cdr u or not eqcar(cadr u,'equal) then
          rederr("Chebyshev_eval: point missing");
       e:=cdr listeval(e,t);
       w:=numr simp caddr cadr u;
       r:= prepf chebeval(e,x,car y,cadr y,w);
    >>
   else if mode = 'int or mode = 'df then
    <<
       e:=cdr listeval(e,t);
       r:= if mode ='int then
                chebint(e,x,car y,cadr y)
           else chebdf (e,x,car y,cadr y);
       r:='list . for each q in r collect prepf q;
    >>;
    switch!-mode!-rd oldmode;
    return r;
  end;

symbolic procedure chebcoeffs(func,x,a,b,n);
  begin
        integer k,j,n1;
	scalar fac,bpa,bma,f,!1pi,c,y,su,nn,half,w;

    x:={x};
    !1pi := rdpi!*();
    nn := dm!:(1/n);
    n1 := n-1;
    half := dm!:(1/2);

     dm!: <<bma:=(b-a)/2; bpa:=(b+a)/2>>;
     w := dm!:(!1pi*nn);
     f:=for k:=0:n1 collect
     <<y := dm!: (rdcos!*(w*(k+half))*bma+bpa);
           evaluate(func,x,{y})
     >>;
        
     dm!: <<fac:=2*nn>>;
     c:=for j:=0:n1 collect
     <<
	w:= dm!:(!1pi*j*nn);
        su:=nil;
        for k:=0:n1 do
          dm!:(su := su+nth(f,iadd1 k) *rdcos!*(w*(k+half)));
        dm!: (fac*su)
     >>;
    return c;
  end;

symbolic procedure chebpol(c,x,a,b);
  begin scalar d,dd,sv,fac,cnst;
        integer n;
    n:=length c;
    d:=for i:=1:n collect nil;
    dd:=for i:=1:n collect nil;
    car dd := nth(c,n);
    for j:=(n-2) step -1 until 1 do
    <<
      for k:=(n-j) step -1 until 1 do
       << sv := nth(d,k+1);
          nth(d,k+1) := dm!:(2*nth(d,k) - nth(dd,add1 k));
          nth(dd,k+1) := sv
       >>;
      sv:=car d;
      car d := dm!:(- car dd + nth(c,add1 j));
      car dd := sv;
    >>;
    for j:=(n-1) step -1 until 1 do
       nth(d,j+1) := dm!:(nth(d,j) - nth(dd,add1 j));
    car d := dm!:(-car dd + car c/2);

    cnst:=dm!:(2/(b-a));
    fac :=cnst;
    for j:=1:(n-1) do
     <<nth(d,add1 j) := dm!:(nth(d,add1 j) * fac);
       fac := dm!:(fac*cnst);
     >>;
    cnst:=dm!:((a+b)/2);
    for j:=0:(n-2) do
      for k:=(n-2) step -1 until j do
        nth(d,add1 k) := dm!:(nth(d,add1 k) - cnst*nth(d,add1 add1 k));

    return reval ('plus. for i:=1:n collect
      {'times,nth(d,i),{'expt,x,i-1}});
  end;

symbolic procedure chebeval(c,xx,a,b,x);
  begin scalar d,dd,y,y2,sv;
        integer m;
    c:=for each q in c collect numr simp q;
    m:=length c;
    y2 := dm!:(2*(y:=(2*x-a-b)/(b-a)));
    for j:=(m-1) step -1 until 1 do
    << sv:=d;
       d:= dm!:(y2*d-dd+nth(c,add1 j));
       dd := sv;
    >>;
    return dm!:(y*d-dd + car c/2);
  end;

symbolic procedure chebint(c,xx,a,b);
  begin scalar su,fac,con,cint,w;
        integer n,jj;
    n := length c;
    c:=for each q in c collect numr simp q;
    fac := 1;
    con := dm!:((b-a)/4);
    cint := for j:=1:(n-2) collect
    << jj:=j+2;
       dm!:(w:= con*(nth(c,j)-nth(c,jj))/j );
       dm!:(su := su + fac*w);
       dm!:(fac := -fac);
       w
    >>;
    cint := append(cint,{dm!: (con*nth(c,sub1 n)/(n-1))});
    dm!:( su := su+fac*nth(c,n));
    cint := (dm!:(su+su)) . cint;
    return cint;
  end;

symbolic procedure chebdf(c,xx,a,b);
  begin scalar con,cder,w;
        integer n,jj;
    n := length c;
    c:=for each q in c collect numr simp q;
    cder := for i:=1:n collect nil;

    nth(cder,n-1) := dm!:(2*(n-1)*nth(c,n));
    for j:=(n-3) step -1 until 0 do
    <<jj:=j+3;
      nth(cder,j+1):=
       dm!:(nth(cder,jj)+2*(j+1)*nth(c,add1 add1 j));
    >>;

    con := dm!:(2/(b-a));

    return for each q in cder collect dm!:(con * q);
  end;

endmodule;


module rungeku;  % Numeric solution of ODE's with Runge-Kutta.
   % Version 2: supporting ODE systems

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1993, all rights resrved

fluid '(!*noequiv accuracy!*);
global '(iterations!* !*trnumeric);

put ('num_odesolve,'psopfn,'rungekuttaeval);

symbolic procedure rungekuttaeval u;
     % interface function;
  begin scalar e,f,x,y,sx,sy,en,d,nd,v,w,q;
        integer n;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,20,6);

     % equations
    e :=car u; u :=cdr u;
    e := if eqcar(e,'list) then cdr e else {e};

     % dependent variable
    q :=car u; u :=cdr u;
    q := if eqcar(q,'list) then cdr q else {q};
    for each yy in q do
    <<if not eqcar(yy,'equal) or not idp cadr yy
     then typerr(yy,"expression `dep. variable=starting value'");
     sy:=caddr yy.sy; y := cadr yy.y;
    >>;
    sy:=reversip sy; y:=reversip y;
    if length sy neq length e then
      rederr "number of equations and variables don't correspond";

      % independent variable
    x :=car u; u :=cdr u;
    if not eqcar(x,'equal) or not idp cadr x
       or null (w:=revalnuminterval(caddr x,t))
        then typerr(x,"expression `indep. variable=interval'");
    sx:=car w; en:=cadr w; x := cadr x;

      % convert expressions to explicit ODE system.
    q := for each yy in y collect {'df,yy,x};
    v:=cdr solveeval list('list. e, 'list . q);
    if cdr v then
      ruksyserr(nil,"cannot convert to explicit ODE system");
    f := for each d in q collect rukufindrhs(d,v);
    return rungekutta1(f,x,y,sx,en,sy);
end;

symbolic procedure rukufindrhs(d,e);
  % Find the righthand side for d in list e.
  if atom e then
      ruksyserr(d," cannot be moved to lhs of system") else
  if eqcar(car e,'equal) and cadr car e = d then
           reval caddr car e else
  if pairp e and eqcar(car e,'list) then  rukufindrhs(d,cdar e)
  else rukufindrhs(d,cdr e);

symbolic procedure ruksyserr(u,m);
  <<if u then writepri(mkquote u,'first);
    writepri(mkquote m,'last);
    rederr "Runge-Kutta failed";
  >>;

symbolic procedure rungekutta1(f,x,y,sx,ex,sy);
    %preparations for rungekutta iteration.
 (begin scalar acc,r,oldmode,!*noequiv;
    integer prec;
    if not memq(dmode!*,'(!:rd!: !:cr))then
       <<oldmode:=t; setdmode('rounded,t)>>;
    !*noequiv:=t; prec:=precision 0;
    sx := force!-to!-dm numr simp sx;
    ex := force!-to!-dm numr simp ex;
    sy := for each z in sy collect force!-to!-dm numr simp z;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    if !*trnumeric then prin2t "starting rungekutta iteration";
    r := rungekutta2(f,x,y,sx,ex,sy,acc);
    if oldmode then setdmode('rounded,nil);
    precision prec;
    if null r then rederr "no solution found";
    return 'list.r;
  end) where !*roundbf=!*roundbf;

symbolic procedure rungekutta2(f,xx,yy,xs,xe,ys,acc);
  % Algorithm for numeric ODE solution
  % f(xx,yy): rhs;
  % x:     independent variable;                    
  % y:     dependent variable; 
  % s:     starting point;
  % e:     endpoint;                                
  % acc:   required accuracy
  % f,yy,ys are vectors (lists).
 dm!:
  begin scalar n,h,h2,h4,d1,d2,d3,d4,x,y,r,w1,w2,vars,dir;
    integer count,st;
    vars := xx.yy;
    x:=xs; y:=ys;
    h:=(xe - xs)/iterations!*;
     dir := h>0; st:= iterations!*;
     % compute initial stepsize
  adapt:
    h2:=h/2; h4:=h/4;
     % full stepsize.
    w1:=rungekuttastep(f,vars,x,y,h,h2);
     % same with half stepsize.
    w2:=rungekuttastep(f,vars,x,y,h2,h4);
    w2:=rungekuttastep(f,vars,x+h2,w2,h2,h4);
      % abs(w2 - w1)
    if normlist list!-list(w2,w1) > acc then
    <<h:=h2; st:=st+st; goto adapt>>;
    if !*trnumeric and not(st=iterations!*) then
      <<prin2
         "*** RungeKutta increased number of steps to ";
        prin2t st>>;
  loop:
    if (dir and x>xe) or (not dir and x<xe) then goto finis;
    r:=('list.x.y) . r;
    count:=add1 count;
    y:= rungekuttastep(f,vars,x,y,h,h2);
      % Advance solution.
    x:=x+h;
    goto loop;
  finis:
    r:= ('list.xx.yy).('list.xs.ys). rungekuttares(reversip r,st);
    return r;
  end;

symbolic procedure rungekuttares(l,st);
   % eliminate intermediate points.
     if st=iterations!* then l else
     << for i:=1:iterations!* collect
         <<for j:=1:m do l:= cdr l; car l>>
     >> where m=st/iterations!*;

symbolic procedure rungekuttastep(f,vars,x,y,h,h2);
  dm!:
   begin scalar d1,d2,d3,d4,q;
             %  f(x,y)
    d1:= list!-evaluate(f,vars,x.y);
             %  f(x+h2,y+h2*d1) 
    d2:= list!-evaluate(f,vars,
          (x+h2). list!+list(y,scal!*list(h2,d1)));
             %  f(x+h2,y+h2*d2)
    d3:= list!-evaluate(f,vars,
          (x+h2). list!+list(y,scal!*list(h2,d2)));
             %  f(x+h,y+h*d3)
    d4:= list!-evaluate(f,vars,
          (x+h).list!+list(y,scal!*list(h,d3)));
      % find y(x+h) = y+h*(d1 + 2*d2 + 2*d3 + d4)/6.
    y:=list!+list(y,
       scal!*list(!:!:quotient(h,6),
         list!+list(d1,
          list!+list(scal!*list(2,d2),
           list!+list(scal!*list(2,d3),
            d4)))));
    return y;
  end;

endmodule;


end;
