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!* boundsvars!*);
symbolic procedure boundsevalrd u;
begin scalar dmode!*;
setdmode('rounded,t);
return boundseval u;
end;
symbolic procedure boundseval u;
(begin scalar db,y,l,f, boundsvars!*;
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);
boundsvars!*:=append(boundsvars!*,{cadr x});
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.
<<boundsvars!*:=for each x in db collect car x;
boundseval1(f,db)>>;
symbolic procedure boundseval1(u,boundsdb!*);
<<if u member boundsknown!* then !*boundsfullcheck := nil else
<< !*boundsfullcheck := t;
if dmode!* = '!:rd!: then 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;
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;
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,n,m;
n := '(nil . 1);
if sqgreaterp(n,ex) then
return boundsquotient('(1 . 1),boundsexpt!-const(b,negsq ex));
if denr ex = 1 and
(fixp(m:=numr ex)
or eqcar(m,'!:rd!:)
and null addf(numr ex,negf(m := reval{'round,m})))
then return boundsexpt!-int(b,m);
if sqgreaterp(n,car b) then
boundserr(nil,"fractional exponent with negative base");
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,n;
n:='(nil . 1);
if sqgreaterp(n,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 boundszero 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;
symbolic procedure boundszero u;
% test: 0 in interval u.
not(minusf numr car u = minusf numr cdr u) or
boundszero1 numr car u or boundszero1 numr cdr u;
symbolic procedure boundszero1 f;
% test standard form f for zero.
null f or zerop f or pairp f and apply(get(car f,'zerop),{f});
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,!1pi,!2pi,!3pi,l,m;
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,prepsq v)) 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,x,ppi;
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,lo,hi,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 member(e, boundsvars!*) 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;
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;
end;