module realroot; % Routines for finding real roots of polynomials,
% using Sturm series, together with iteration.
% Author: Stanley L. Kameny <stan_kameny@rand.org>.
% Version and Date: Mod 1.96, 30 March 1995.
% Copyright (c) 1988,1989,1990,1991,1992,1993,1994,1995.
% Stanley L. Kameny. All Rights Reserved.
Comment modules bfauxil, bfdoer, bfdoer2, complxp, allroot and rootaux
needed also;
exports accupr1, bfnewton, isolatep, schinf, schplus, sgn1, sturm,
sturm0, uniroots;
imports !!mfefix, abs!:, accupr, accuroot, allroots, automod, bdstest,
bfabs, bfdivide, bfeqp, bfleqp, bfloat, bfloatem, bfmax,
bfminus, bfminusp, bfplus, bfrlmult, bfsgn, bfsqrt, bfzp,
ceillog, ckpzro, cpxp, csep, difbf, divbf,
domainp, dsply, eqcar, equal!:, errach, geq, getprec, gfdiff,
gffinitr, gfgetmin, gfrl, gfrootfind, gfsqfrf, gfstorval,
greaterp!:, lastpair, leq, lprim, minbnd1, minprec, mk!*sq,
multroot,
neq, nwterr, nwterrfx, outecho, pconstr, plubf, powerchk,
r2bf, r2flbf, ratdif, ratleqp, ratlessp, ratmax, ratmean,
ratmin, ratminus, ratplus, realrat, rerror, rl2gf, rlval,
round!:mt, sch, schnok, setprec, sgn, stuffr, sturm1, timbf,
trmsg1, trmsg10, trmsg2, trmsg3, trmsg4, trmsg6, trmsg7,
trmsg8, xclp;
global '(!!nfpd !!flim bfhalf!* max!-acc!-incr bfone!* rlval!#);
fluid '(!*gfp !*xnlist !*intp tht!# !*strm lims!# mltr!# pfactor!#
prec!# !*rvar acc!# !*xo !1rp accm!# !*xn intv!# sh!# rprec!# rlrt!#
prm!# pfl!# acfl!# pgcd!#);
fluid '(!*trroot !*bftag !*compxroots !*msg);
flag('(positive negative infinity),'reserved);
global '(limlst!# lm!#);
limlst!# := '(positive negative infinity (minus infinity));
symbolic procedure isolatep p;
begin scalar n,q,zr,a,b,c,ril,va,vb,vc,v0,w,elem,l,u,i,j,lr,ur,
xcli,xclj,ol,ou;
if null sturm p or schinf(-1)-schinf 1=0 then go to ret;
% limits +/-1.0001*maxbound p to give working room for rootfind.
n := car(q := car !*strm);
l := ratminus(u := realrat bfrlmult(1.0001,bfmax p));
if (zr := l=u) and (lr := l) and not lims!# then go to zrt;
if lims!# then
<<i := car lims!#; if cdr lims!# then
<<j := cadr lims!#; % both limits given.
if i eq 'minfty then xcli := t else
<<if xclp i then <<xcli := t; i := cdr i>>;
l := ratmax(l,i)>>;
if j eq 'infty then xclj := t else
<<if xclp j then <<xclj := t; j := cdr j>>;
u := ratmin(u,j)>>;
if zr then if ratlessp(u,l) then go to ret
else go to zrt;
if sgn1(q,l)=0 then % root at l.
<<ol := offsetr(n,l);
if xcli then l := ratplus(l,ol) else lr := l>>;
if l neq u and sgn1(q,u)=0 then % root at u.
<<ou := offsetr(n,u);
if xclj then u := ratdif(u,ou) else ur := u>> >>
else if zr then go to ret else
<<if sgn1(q,ol := realrat 0)=0 then ol := offsetr(n,ol);
if i<0 then u := ratminus ol % negative roots.
else l := ol>> >>; % positive roots.;
n := (va := sch l+if lr then 1 else 0)-(vb := sch u);
trmsg4(n);
if n=0 then go to ret;
if n=1 then ril := list list(l,u)
else for j:=1:n do ril := nil . ril;
v0 := vb+n-1;
if lr then
<<stuffrt(l,u,lr,ol,v0,va,vb,nil,ril);
l := ratplus(lr,ol); va := va-1>>;
if ur then
<<stuffrt(l,u,ur,ou,v0,va,vb,nil,ril);
u := ratdif(ur,ou); vb := vb+1>>;
w := list list(l,u,va,vb);
if n>1 then while w do
<<elem := car w; w := cdr w; a := car elem; b := cadr elem;
va := caddr elem; vb := cadddr elem; c := ratmean(a,b);
if sgn1(q,c)=0 then % c is a root.
w := stuffrt(a,b,c,offsetr(n,c),v0,va,vb,w,ril) else
% root not found; stuff isolating interval and update work.
<<vc := sch c;
if va = vc+1 then <<stuffr(v0-vc,list(a,c),ril)>>;
if va > vc+1 then w := list(a,c,va,vc) . w;
if vb = vc-1 then <<stuffr(v0-vb,list(c,b),ril)>>;
if vb < vc-1 then w := list(c,b,vc,vb) . w>> >>;
ril := for each i in ril collect (car i) . cadr i;
ret: return ril;
zrt: return list (lr . lr) end;
symbolic procedure stuffrt(a,b,c,m,v0,va,vb,w,ril);
begin scalar vcm,vcp; % stuff root and update work.
vcm := 1+(vcp := sch ratplus(c,m));
stuffr(v0-vcp,list(c,c),ril);
if va = vcm+1 then stuffr(v0-vcm,list(a,ratdif(c,m)),ril);
if va > vcm+1 then w := list(a,ratdif(c,m),va,vcm) . w;
if vb = vcp-1 then stuffr(v0-vb,list(ratplus(c,m),b),ril);
if vb < vcp-1 then w := list(ratplus(c,m),b,vcp,vb) . w;
return w end;
symbolic procedure offsetr(n,r);
realrat if n=1 then 1 else minbnd1(!*gfp,mk!*sq r);
symbolic procedure sturm p;
<<if cpxp (p := gffinitr p) then
<<p := car csep p; if not atom p then p := bfloatem p>>;
if not atom p then sturm1(!*gfp := p)>>;
put('sturm,'psopfn,'sturm0);
symbolic procedure sturm0 p;
<<p := sturm ckprec car p; restorefl();
'list . for each a in p collect if atom a then a else 'list . a>>;
symbolic procedure sgn1(p,r); if atom p then sgn p else
% Evaluate sign of one sturm polynomial for rational r=(u . d)
begin scalar m,c,u,d; u := car r; d := cdr r;
c := 0; m := 1; p := cdr p;
repeat <<c := m*car p + u*c; m := m*d>> until null(p := cdr p);
return sgn c end;
symbolic procedure r2flimbf x;
if acc!#<=!!flim then r2flbf x else r2bf x;
symbolic procedure rootfind(p,i);
% finds real roots in either float or bigfloat modes.
% p is in gfform. i is a pointer to a rational interval pair;
begin scalar p1,p2,px,x1,x2,x3,x0,nx,xr,fg,n,s,sh;
scalar xd,xe,qt,xnlist,pf,pf0; integer m,tht!#;
n := caar lastpair p; !*xnlist := nil;
if car i=cdr i then
<<nx := r2flbf cdr i; go to lg4>>;
xr := ratmean(car i,cdr i);
if !*trroot then
<<write "real root ",list(r2flimbf car i,r2flimbf cdr i);
terpri()>>; trmsg8();
if ratlessp(cdr i,car i) then
errach "lx > hx ***error in roots package";
movebds(i,xr,sh!# := sh := sgn1(!*intp,cdr i));
p2 := gfdiff(p1 := gfdiff p);
lag0: if bndtst (px := rlval(p,nx := r2flbf xr)) then go to tht
else if bfzp px then go to lg4;
lag: % check for proper slope at nx.
if bndtst (x1 := rlval(p1,nx))
or (s := bfsgn x1) neq sh then go to tht;
% if lag not converging, go to newt.
pf := bfabs px; if pf0 and bfleqp(pf0,pf) then go to newt;
gfstorval(pf,nx); x1 := bfabs x1;
if bndtst (x3 := rlval(p2,nx)) then go to tht;
% bigfloat computations: is newton cheaper?
if fg and
<<qt := divbf(px,x1);
xe := timbf(qt,timbf(qt,(divbf(x3,x1))));
equal!:(nx,plubf(nx,timbf(bfhalf!*,xe)))>>
then go to newt;
% check whether laguerre iteration will work.
x2 := difbf(bfrlmult(n-1.0,timbf(x1,x1)),
bfrlmult(n,timbf(px,x3)));
if bfminusp x2 then go to tht;
% nx has met all tests, so continue.
x0 := nx;
xd := divbf(bfrlmult(-n*s,px),
plubf(x1,bfsqrt(bfrlmult(n-1,x2))));
nx := plubf(x0,xd);
lg3: fg := t;
if ratlessp(xr := realrat nx,car i) or ratlessp(cdr i,xr)
then go to tht;
if bndtst (px := rlval(p,nx)) then go to tht;
movebds(i,xr,sh); trmsg2 ('lag,nx,px);
if bdstest i then go to ret;
if bfzp px then go to lg4;
if bfeqp(nx,x0) then <<trmsg3('lag,nx); go to ret>>;
if xnlist and member(nx,xnlist) then go to newt;
xnlist := nx . xnlist; pf0 := pf;
if(m := m+1)<10 or
<<m := 0;
equal!:(bfone!*,round!:mt(divbf(bfloat nx,bfloat x0),13))>>
then go to lag;
tht: nx := tighten(i,p,pf,sh); m := 0;
if !*xnlist then
<<pf0 := nil;
movebds(i,xr := ratmean(car i,cdr i),sh); go to lag0>>;
lg4: trmsg1('lag,nx);
ret: !*xnlist := nil; if not nx then trmsg10 'lag; go to ret2;
newt: nx := bfnewton(p,p1,gfgetmin(),i,4);
ret2: !*xn := rl2gf nx; return nx end;
global '(tentothetenth!*!*);
tentothetenth!*!* := normbf i2bf!: 10000000000;
symbolic procedure bndtst x;
greaterp!: (abs!: x, tentothetenth!*!*);
symbolic procedure movebds(i,xr,sh);
if sgn1(!*intp,xr)=sh then rplacd(i,xr) else rplaca(i,xr);
symbolic procedure tighten(i,p,pf,sh);
begin scalar j,x0,nx,px,sn,x;
nx := car i;
tht0: j := 4;
tht1: x0 := nx; nx := ratmean(car i,cdr i);
if (sn := sgn1(!*intp,nx))=0 then
<<x := r2flbf nx;trmsg1 ('tht,x); go to ret>>;
if 0=car ratdif(nx,x0) then
<<x := r2flbf nx;trmsg3 ('tht,x); go to ret>>;
if sn=sh then rplacd(i,nx) else rplaca(i,nx);
if (sn := bdstest i) then <<x := r2flbf sn; go to ret>>;
if (j := j-1)>0 then go to tht1;
if bndtst (px := rlval(p,x := r2flbf nx)) then
<<j := 4; go to tht1>>;
gfstorval(bfabs px,x);
trmsg2('tht,x,px);
if bfzp px then go to ret
else if pf and bfleqp(pf,bfabs px) then go to tht0
else return x;
ret: !*xnlist := nil; return x end;
symbolic procedure rtsreal(p,s);
% Finds real roots of univariate square-free real polynomial p, using
% sturm series, isolater and rootfind.
begin scalar acr,acs,n,q,r,x,y,!*strm,pr,apr,!*bftag,pfl!#,
acfl!#,xout,x1;
integer accn,accn1,accm!#,prec!#,prm!#; pr := getprec();
!*bftag := rlrt!# := t; pgcd!# := not s;
r := isolatep p; % r is a list of rational number pairs.
if null r then go to ret;
if (n := caar lastpair p)>1 then go to gr1;
y := rootrnd gfrl gfrootfind(p,nil);
if pfactor!# then
<<y := accupr1(y,p); y := (rootrnd car y) . cdr y>>;
% note that rlval!# was set by the last operation of rootrnd.
xout := {if s then (mkdn rlval!#) . acc!#
else if pfactor!# then y
else y . acc!# % this can't happen
};
if !*trroot then terpri(); go to ret;
gr1: !*xo := rl2gf 0;
q := r; acs := acc!#;
lag: % increase accuracy for this root and the next root if current
% accuracy is not sufficient for the interroot interval.
if cdr q then % no test if this is the last real root.
<<setprec acs;
while schnok q do setprec (getprec()+1);
accn1 := getprec()>>;
acc!# := max(acs,accn,accn1);
accn := if accn1>acs then accn1 else 0;
setprec max(rprec!#,acc!#+2);
y := rootfind(p,intv!# := car q); apr := t;
if null y then rerror(roots,8,"Realroots abort");
acc: y := accuroot(gfrl !*xn,p,!*xo);
% if acc!# is insufficient for this root, for any reason,
% increase accuracy and tighten.
if apr then
<<if (acr := accupr(p,!1rp,!*xn))>acc!# then acc!# := acr
else if acr<=acc!# then <<acc!# := acr; apr := nil>>;
go to acc>>;
xout := ((x1 := if s then mkdn rlval!# else y) . acc!#) . xout;
% x is root list. Check for equal roots should fail!
if x and x1=car x then rooterr x1;
x := x1 . x;
dsply y;
acc!# := acs;
if (q := cdr q) then <<accn1 := 0; go to lag>>;
ret: setprec pr; return reverse xout end;
symbolic procedure lval x; if xclp x then cdr x else x;
symbolic procedure lpwr(l,m);
if eqcar(l,'list) then 'list . lpwr(cdr l,m)
else if atom l then l else ((car l)**m) . ((cdr l)**m);
symbolic procedure schnok r;
%true if precision is inadequate to separate two adjacent real roots.
(l neq h and (sch l neq sch r2flbf2r l or sch h neq sch r2flbf2r h))
where l=caar r,h=cdar r;
symbolic procedure limchk x;
<<!!mfefix();
if null (x := for each y in x collect
if member(y,limlst!#) then y
else if eqcar(y := reval y,'list)
then 'list . list limchk1 cadr y
else limchk1 y) then nil
else if x and not cdr x then
if car x eq 'positive then list 1
else if car x eq 'negative then list (-1) else limerr()
else <<x := mkratl x; limchk2(car x,cadr x)>>>>;
symbolic procedure limchk1 y;
if errorp(y := errorset!*({'a2rat,mkquote y},nil))
then rerror(roots,5,"Real root function limits must be real")
else car y;
symbolic procedure limchk2(a,b);
<<if member(a,l) and member(b,l) then if a neq b then nil else limerr()
else if member(a,limlst!#) then
if member(b,limlst!#) then limerr() else limchk2(b,a)
else if member(b,limlst!#) then
if b eq 'negative then list('minfty,mkxcl a)
else if b eq 'positive then list(mkxcl a,'infty)
else if b eq 'infinity then list(a,'infty) else list('minfty,a)
else if ratv b=ratv a and (xclp a or xclp b) then t
else if ratlessp(ratv b,ratv a) then list(b,a) else list(a,b)>>
where l = cddr limlst!#;
symbolic procedure limerr;
rerror(roots,6,"Illegal region specification");
symbolic procedure ratv a; if xclp a then cdr a else a;
symbolic procedure a2rat x;
if numberp x then x . 1
else if atom x then limerr()
else if eqcar(x,'quotient) then
((if numberp n then n
else if eqcar(n,'minus) then - cadr n
else rerror(roots,10,"illegal limit")) where n=cadr x) . caddr x
else if car x eq '!:rn!: then cdr x
else
((if car x memq domainlist!* and y then cdr(apply1(y,x))
else limerr())
where y=get(car x,'!:rn!:));
symbolic procedure rlrootno a;
<<mltr!# := t; lims!# := limchk cdr a; a := ckprec car a;
a := rlrtno2 if lims!#=t then 0 else a;
restorefl(); mltr!# := nil; a>>;
put('rlrootno,'psopfn,'rlrootno);
symbolic procedure realroots a;
<<lims!# := limchk cdr a;
uniroots(if lims!#=t then 0 else car a,0)>>;
put('realroots,'psopfn,'realroots);
symbolic procedure isolater p;
<<mltr!# := t; lims!# := limchk cdr p; p := ckprec car p;
p := isolatep if lims!#=t then 0 else p;
restorefl(); mltr!# := nil; outril p>>;
put('isolater,'psopfn,'isolater);
symbolic procedure mkratl l; for each a in l collect
if member(a,limlst!#) then a else
if eqcar(a,'list) then
if member(a := cadr a,limlst!#) then a else mkxcl a
else a;
symbolic procedure exclude x; {'list, x};
symbolic operator exclude;
endmodule;
end;