% ----------------------------------------------------------------------
% $Id: ofsfanuex.red,v 1.9 2003/02/02 22:13:12 seidl Exp $
% ----------------------------------------------------------------------
% Copyright (c) 2001 A. Dolzmann, A. Seidl, and T. Sturm
% ----------------------------------------------------------------------
% $Log: ofsfanuex.red,v $
% Revision 1.9 2003/02/02 22:13:12 seidl
% Bug in verbose output fixed.
%
% Revision 1.8 2003/01/29 17:35:28 seidl
% Even second-level verbose output now depends on switch rlverbose.
%
% Revision 1.7 2002/05/28 13:05:04 seidla
% constp added.
%
% Revision 1.6 2002/02/19 13:36:48 seidla
% Minor changes.
%
% Revision 1.5 2002/01/25 15:30:43 seidla
% Further optimization in ratpoly_subrat and aex_remseq. New function
% aex_pp for primitive part.
%
% Revision 1.4 2002/01/23 18:29:26 seidla
% Optimization in ratpoly_subrat and aex_remseq.
%
% Revision 1.3 2002/01/16 16:10:09 seidla
% Small bug resolved.
%
% Revision 1.2 2002/01/16 15:00:37 sturm
% Changed module name.
%
% Revision 1.1 2002/01/16 13:03:48 seidla
% Imported CAD from rlprojects.
%
% Revision 1.31 2002/01/16 11:23:22 seidla
% Sorting of sturm chains for root isolation enabled. Optimization in
% aex_quotrem. Optimization in aex_rem turned out not to be a good idea.
%
% Fine tuning in many places.
%
% calculation of sign changes of sturm chains now at infinity possible.
% this is used now to check if all roots are found. gcd does not
% normalize. reduce added for qqi in aex_quotrem.
%
% New version of anuex.
%
% Revision 1.30 2001/12/11 17:43:22 seidla
% minimization in aex_red. as reimplementation of anuex has started,
% support of this branch will cease soon.
%
% Revision 1.29 2001/11/26 16:50:22 seidla
% now it compiles.
%
% Revision 1.28 2001/08/20 09:55:30 seidla
% Final version for diploma work.
%
% Revision 1.27 2001/08/12 18:40:39 seidla
% Removed new bug in aex_pseudorem1. Changes to aes_remseq, aex_sqfree,
% aex_sgn, aex_nullp, aex_findroots wrt base cases. Now some examples
% run considerably faster. Made some comments more precise.
%
% Revision 1.26 2001/08/11 19:47:06 sturm
% Working. Changed a lot. Currently at aex_isoroot(ae,m,r,sc);
%
% Revision 1.25 2001/08/11 17:20:41 seidla
% Various little changes.
%
% Revision 1.24 2001/08/10 10:14:46 seidla
% Cleaned up a lot and moved unused functions to anuexplus.red.
%
% Revision 1.23 2001/08/07 21:23:57 sturm
% Changed some minor things.
%
% Revision 1.22 2001/08/07 12:56:24 seidla
% aex_mapmkaex0d and anu_findrootsoflist from package OFSFCAD moved to
% here. anu_evalsignf moved to package OFSFCAD as acad_evalsignf.
%
% Revision 1.21 2001/06/14 21:17:03 sturm
% Introduced fluids anuex_rcsid!* and anuex_copyright!*.
% Made a hack in anu_evalsignf to overcome the non-implemented case where
% some cylinder consists of only one cell. (The test point is then still
% "nil").
%
% Revision 1.20 2001/06/14 19:36:53 seidla
% New function anu_evalsignf.
%
% Revision 1.19 2001/06/13 11:37:21 seidla
% aex_findroots rewritten. new: aex_pseudoqrem, aex_gcd, aex_sqfree.
% aex_pseudoqrem tested with aex_pseudoqremtest. speed-up for aex_comp
% by using ssqarefree part of gcd instead of multiplication. new switch
% anuexdebug. anu_mergesort sorts out representations of the same number.
%
% Revision 1.18 2001/06/12 13:43:59 sturm
% Added a variant of anu_compdifferent under the name anu_compdifferent!-ts.
%
% Revision 1.17 2001/05/30 14:58:29 seidla
% Sorting and comparison of algebraic numbers was rewritten. We refer to
% the example cox14v. For algebraic numbers which are known to be
% different, the base phase is now faster by factor 9. For algebraic
% numbers which are not necessarily different, sorting is now feasible,
% i.e. the computing time was reduced from >17h to 1000s. The changes
% made together with the time (elsie=ultra1, bender=ultra10):
% % Time: 41600 ms plus GC time: 2300 ms on elsie *** starting point
% % Time: 37630 ms plus GC time: 2010 ms on elsie *** unnecessary item
% removed in anu_refine1
% % Time: 32090 ms plus GC time: 1770 ms on elsie *** optimization in
% aex_remseq and anu_comp
% % Time: 20220 ms plus GC time: 1110 ms on elsie *** my new mergesort
% % Time: 17990 ms plus GC time: 1100 ms on elsie *** sturmchain handed
% down in anu_compdifferent
% % Time: 17860 ms plus GC time: 1080 ms on elsie *** now everything runs
% over anu_comp
% % Time: 11330 ms plus GC time: 640 ms on elsie *** anu_comp calls
% anu_compdifferent with smaller intervals
% % Time: 11210 ms plus GC time: 630 ms on elsie *** apply reversip
% in anu_mergesort1 to avoid trailing small lists
% % Time: 4696600 ms plus GC time: 4990 ms on bender *** now the base
% phase ran the first time with anuexdifferentroots off.
% % Time: 7900 ms plus GC time: 430 ms on elsie *** in anu_comp intervals
% are changed.
% % Time: 4630 ms plus GC time: 230 ms on elsie *** anu_refineip is used
% % Time: 1999180 ms plus GC time: 3210 ms on bender ****
% anuexdifferentroots off
% anu_check checks validy of algebraic numbers. aex_findroots is still
% dodgy and has to be rewritten. After mergesort, still not all
% intervals are distinct due to some remaining cases in anu_comp. new
% switches anuexdifferentroots (to switch between sorting of different
% and not necessarily different numbers) and anuexmergelookahead.
%
% Revision 1.16 2001/05/25 08:34:08 seidla
% Sorting now uses anu_compdifferent. This works out for the basis phase,
% but it is too slow. Hence sorting has to be rewritten.
%
% Revision 1.15 2001/05/24 15:24:53 sturm
% Alternative implementation for anu_mergesort.
% Removed begin inside progn to the toplevel in aex_sgn.
% BUG: infinite recursion in aex_sgn.
%
% Revision 1.14 2001/05/24 13:55:16 seidla
% nu_mergesort now returns a list of algebraic numbers with pairwise
% distinct isolating interval. new functions anu_distinguish,
% anu_distinguishlist, anu_refine.
%
% Revision 1.13 2001/05/24 12:07:09 sturm
% bug in aex_findroots fixed.
%
% Revision 1.12 2001/05/24 09:09:34 seidla
% changes to aex_rem. aex_nullp and aex_sgn rewritten.
%
% Revision 1.11 2001/05/24 08:51:16 sturm
% Corrected specification comment for <CONTEXT>.
%
% Revision 1.10 2001/03/29 16:47:11 seidla
% New procedure anu_fromratnex.
%
% Revision 1.9 2001/03/27 12:11:14 seidla
% New procedures anu_less and anu_leq. Code cosmetics: No capital
% letters in procedure variables anymore. New procedure
% anuex_initialize.
%
% Revision 1.8 2001/03/26 21:00:45 seidla
% aex_cauchybound works now for all cases. bug in int_add and anu_2aex
% found. Code cosmetics: spaces around commas and capital letters in
% procedure names removed.
%
% Revision 1.7 2001/03/26 13:16:34 sturm
% Exported some more procedures.
%
% Revision 1.6 2001/03/26 12:04:04 seidla
% Exported functions added. anu_mkAfromQ renamed to anu_mkfromaex,
% aex_mkQ0d renamed to aex_mkaex0d, aex_getL renamed to aex_getctx and
% anu_getL renamed to anu_getctx.
%
% Revision 1.5 2001/03/22 15:47:31 sturm
% Removed TeX in first sentence of comments.
%
% Revision 1.4 2001/03/22 14:09:52 seidla
% Code reformatted.
%
% Revision 1.3 2001/03/21 22:17:18 seidla
% Code cosmetics: all comments have been rewritten to a common style.
% A part of the code was reformatted, but not all.
% Two bugs were found in aex_findrootsiniv.
%
% Revision 1.2 2001/03/21 12:34:37 sturm
% Did some first code formatting.
%
% Revision 1.1 2001/03/15 12:30:34 seidla
% Initial check-in.
%
% ----------------------------------------------------------------------
% file: anuex.red
% author: andreas seidl
%
% date: 13.03.2001
% sortierung fertiggestellt. sicherheitskopie gemacht. ran.red in anuex.red
% umbenannt. sublist in context umbenannt.
% date 12.03.2001
% findrootsiniv fertiggestellt. fehler im testfile: q3. sortierung angefangen.
% date 09.03.2001
% findroots begonnen
% ein vermeindlicher bug in aex_sgn stellte sich als fehler im testfile
% heraus.
% date 08.03.2001
% aex_containment und aex_cauchybound (falls a_m<>0) fertiggestellt.
% findroots begonnen.
% date 26.02.2001
% aex_nullp, aex_sgn. sicherheitskopien von ran.red, anuex.test, others.test.
% date 23.02.2001
% f"ur aex_nullp und aex_sgn bereits rest bei pseudodivision und minimize
% implementiert.
% date 16.02.2001
% rekursive datenstruktur notwendig. konstruktoren und zugriffsfunktionen
% f"ur anu und aex.
% date: 09.02.2001
% sturmchain und darauf aufbauend ran_sgn implementiert
% date: 08.02.2001
% kapitel package.tex begonnen, wo datentypen genauer spezifiziert sind.
% ran_add nach ran_add2 und ran_mult nach ran_mult2 umbenannt
% redundand datatype ran_poly_* gel"oscht.
% _int_ von _num_ auf _rat_ umgestellt, entsprechende "anderungen in
% ratpoly_sub. damit scheint ran_nullp nun zu funktionieren. weitere
% tests sind notwendig.
% date: 03.02.2001
% ran_sqd nach ratpoly umbenannt, quadrat-frei
% date: 29.01.2001
% datentyp sqd
% date: 23.01.2001
% rationals are needed to allow for polynomial division.
% hence ran_poly2sf rewritten.
% date: 22.01.2001
% ran_sfdeg: numberp replaced by domainp
% date: 19.01.2001
% datentyp poly und ran
% date: 09.01.2001
% datentyp int
lisp <<
fluid '(ofsfanuex_rcsid!* ofsfanuex_copyright!*);
ofsfanuex_rcsid!* := "$Id: ofsfanuex.red,v 1.9 2003/02/02 22:13:12 seidl Exp $";
ofsfanuex_copyright!* := "(c) 2001 by A. Dolzmann, A. Seidl, T. Sturm"
>>;
module ofsfanuex;
% ANUEX - a package for primitive or hierarchical representation of
% algebraic numbers and expressions in algebraic numbers (algebraic
% polynomials). the features include incremental real root isolation
% and factorisation.
procedure ofsf_anuexverbosep();
% CAD verbose predicate.
!*rlverbose and !*rlanuexverbose;
procedure type_of(te);
% Type of. [te] is a typed expressions
car te;
% module generic;
procedure generic_sort(al,sortfn);
% Sort. [l] is a list of some type ALPHA, [sortfn] is a
% function that takes two arguments [a1], [a2] of type ALPHA and
% returns 1, if [a1] is bigger than [a2], -1, if [a2] is bigger
% than [a1] and 0 if they are equal. returns a list of type ALPHA.
begin scalar all;
if null al or null cdr al then return al;
% now there are at least 2 elements in al.
all := for each a in al collect {a};
repeat <<
all := generic_doublesort(all,sortfn)
>> until null cdr all;
return car all % the only element of all is the list of sorted roots
end;
procedure generic_doublesort(all,sortfn);
begin scalar res;
% as long as there at least two elements in all:
while not (null all or null cdr all) do <<
res := generic_doublesort1(car all,cadr all,sortfn) . res;
all := cddr all;
>>;
if all then res := car all . res;
return res;
end;
procedure generic_doublesort1(al1,al2,sortfn);
% Algebraic number doublesort. Takes two sorted lists. Returns a
% sorted list.
begin scalar rl,cmp;
if null al1 then
return al2;
if null al2 then
return al1;
while al1 and al2 do <<
cmp := apply2(sortfn,car al1,car al2);
if cmp eq 0 then <<
rl := car al1 . rl; al1 := cdr al1; %al2 := cdr al2
>> else if cmp eq -1 then << % car al1 < car al2
rl := car al1 . rl; al1 := cdr al1
>> else <<
rl := car al2 . rl; al2 := cdr al2
>>
>>;
return nconc(nconc(reversip rl,al1),al2)
end;
% module iv;
% Interval
%DS
% <INT> ::= (<RAT> . <RAT>)
% <RAT> ::= <SQ> over the integers
procedure iv_mk(l,r);
% Interval make. [l] and [r] are RAT's. Returns an INT.
l . r;
procedure iv_mk4(n1,d1,n2,d2);
iv_mk(rat_mk(n1,d1),rat_mk(n2,d2));
procedure iv_lb(i);
% Interval left bound. [i] is an INT. Returns a RAT. Returns the
% left bound of [i].
car i;
procedure iv_rb(iv);
% Interval right bound. [i] is an INT. Returns a RAT. Returns the
% right bound of [i].
cdr iv;
procedure iv_print(iv);
<< prin2 "]"; rat_print iv_lb iv; prin2 ",";
rat_print iv_rb iv; prin2 "[";
>>;
procedure iv_neg(i);
% Interval negate. [i] is an INT. Returns an INT. Returns $-[i]$,
% i.e., the interval mirrored at 0.
iv_mk(negsq iv_rb I,negsq iv_lb I);
procedure iv_add(iv1,iv2);
% Interval addition. [iv1] and [iv2] are INT. Returns an INT.
iv_mk(addsq(iv_lb iv1,iv_lb iv2),addsq(iv_rb iv1,iv_rb iv2));
procedure iv_mapadd(ivl);
% Interval map addition. [ivl] is a list of INT. Returns an interval.
% Recommendation: [ivl] is non-empty
if null ivl then
iv_mk(rat_fromnum 0,rat_fromnum 0)
else
iv_add(car ivl,iv_mapadd cdr ivl);
procedure iv_mult(iv1,iv2);
% Interval multiplication. [iv1] and [iv2] are INT. Returns an INT.
iv_mk(rat_min4(rr,rl,lr,ll),rat_max4(rr,rl,lr,ll)) where
rr = multsq(iv_rb iv1,iv_rb iv2),
rl = multsq(iv_rb iv1,iv_lb iv2),
lr = multsq(iv_lb iv1,iv_rb iv2),
ll = multsq(iv_lb iv1,iv_lb iv2);
procedure iv_tothen(iv,n);
% Interval to the n. [iv] is an interval, [n] a positive natural
% number. Returns an INT.
if n=0 then
iv_mk(rat_mk(1,1),rat_mk(1,1))
else
iv_mult(iv,iv_tothen(iv,n-1));
procedure iv_containszero(iv);
% Interval contains zero. Returns [t] or [nil].
if rat_leq(iv_lb iv,rat_0()) and rat_leq(rat_0(),iv_rb iv) then
t
else
nil;
procedure iv_comp(iv1,iv2);
% Interval compare intervals. [iv1], [iv2] are intervals. Returns -1,
% if [iv1] is below [iv2], 1, if [iv2] is above [iv1] and 0, if the
% intersection of [iv1] and [iv2] is non- empty. [iv1], [iv2] are
% regarded as open.
if rat_leq(iv_rb iv1,iv_lb iv2) then
-1
else if rat_leq(iv_rb iv2,iv_lb iv1) then
1
else 0;
procedure iv_maxabs(iv);
% Interval maximal absolut value. [iv] is an INT. Returns a RAT.
rat_max(rat_abs iv_lb iv,rat_abs iv_rb iv);
procedure iv_minabs(iv);
% Interval minimal absolut value. [iv] is an INT. Returns a RAT.
if iv_containszero iv then
rat_0()
else
rat_min(rat_abs iv_lb iv,rat_abs iv_rb iv);
procedure iv_minus(iv1,iv2);
% Returns a list of IV, resulting from subtracting [iv2] from [iv1].
nconc(
if rat_less(iv_lb iv1,iv_lb iv2) then
{iv_mk(iv_lb iv1,rat_min(iv_lb iv2,iv_rb iv1))},
if rat_less(iv_rb iv2,iv_rb iv1) then
{iv_mk(rat_max(iv_rb iv2,iv_lb iv1),iv_rb iv1)});
procedure iv_minuslist(iv1,ivl2);
% ivl2 is a list of distinct intervals. Returns a list of IV.
if null ivl2 then
{iv1}
else
iv_listminuslist(iv_minus(iv1,car ivl2),cdr ivl2);
procedure iv_listminuslist(ivl1,ivl2);
% ivl1, ivl2 are each lists of distinct intervals.
for each iv1 in ivl1 join iv_minuslist(iv1,ivl2);
% module num;
%DS
% <NUM> ::= integer
procedure num_even(n);
n eq n/2*2;
procedure num_sgnch(l);
% Integer number sign changes. [l] is a list of NUM. Returns a
% NUM.
num_sgnch1(for each n in l join if sgn n = 0 then nil else n . nil);
procedure num_sgnch1(l);
% Integer number sign changes 1. [l] is a list of positive or
% negative NUM. Returns a non-negative NUM. Counts the number
% of sign changes.
if null l or null cdr l then
0 % in a list with at most one element there are no sign changes
else if sgn car l eq sgn cadr l then
num_sgnch1 cdr l
else
1 + num_sgnch1 cdr l;
procedure num_sortfn(n1,n2);
% for use in generic_sort.
sgn(n1-n2);
% module numpoly;
%DS
% <NUMPOLY> ::= <SF> over the integers
procedure numpoly_null;
% Integer number polynomial null. Returns a NUMPOLY.
nil;
procedure numpoly_nullp(p);
% Integer number polynomial null predicate. [p] is a
% NUMPOLY. Returns [t] or [nil]. Double null is allowed for.
null p or p eq 0;
procedure numpoly_lc(p);
% Integer number polynomial leading coefficient. [p] is a
% NUMPOLY. Returns a NUMPOLY. Caveat: extended semantics wrt. lc
% for sf's.
if null p then
p
else if numberp p then
p
else
lc p;
procedure numpoly_red(p);
% Integer number polynomial reductum. p is a numpoly. Returns a
% numpoly.
if null p or numberp p then
numpoly_null()
else
red p;
procedure numpoly_mvartest(p,x);
% Integer number polynomial main variable test. [p] is a NUMPOLY,
% [x] is an identifier. Returns [t] or [nil]. Checks, whether [x]
% is the main variable of [p]
if null p then
nil
else if numberp p then
nil
else
(mvar p eq x);
% REMARK: the following holds: p a numpoly in x. then:
% p == (if mvartest(p,x) then mult(lc(p),xtothen('x,ldeg p))+red p
procedure numpoly_ldeg(p);
% Integer number polynomial leading degree. [p] is a
% NUMPOLY. Returns -1 or a positive integer.
if numpoly_nullp p then
-1
else if numberp p then
0
else
ldeg p;
procedure numpoly_factorize p;
fctrf p;
% module rat;
% Rational numbers.
%DS
% <RAT> ::= <SQ> over integers
% datatype rat: rational numbers, represented as standard quotients of
% integers.
procedure rat_mk(n,d);
% Rational number make. [n] and [d] are NUM. Returns a RAT.
quotsq(simp n,simp d);
% n . d;
procedure rat_numr(q);
numr q;
procedure rat_denr(q);
denr q;
procedure rat_normdenr(q);
<< if !*rlanuexdebug then if minusf rat_denr q then
prin2 "***** rat_normdenr: denr is negative";
rat_numr q ./ 1
>>;
procedure rat_print r;
<< prin2 numr r; prin2 "/"; prin2 denr r; >>;
procedure rat_fromnum(n);
% Rational number from integer number. [n] is a NUM. Returns a RAT.
simp n;
procedure rat_neg(p);
% Rational number negation. [p] is a RAT.
negsq(p);
procedure rat_add(p,q);
% Rational number addition. [p], [q] are RAT's. Returns a RAT.
addsq(p,q);
procedure rat_minus(p,q);
% Rational number minus. [p], [q] are RAT's. Returns a RAT.
addsq(p,negsq q);
procedure rat_mapadd(pl);
% Rational number map addition. [pl] is a list of RAT. Returns a RAT.
if pl=nil then %%%
rat_0()
else
rat_add(car pl,rat_mapadd cdr pl);
procedure rat_mult(p,q);
% Rational number multiplication. [p], [q] are RAT's. Returns a RAT.
multsq(p,q);
procedure rat_quot(p,q);
% Rational number quotient. [p], [q] are RAT's, the latter has to
% be non-zero. Returns a RAT.
quotsq(p,q);
procedure rat_nullp(q);
% Rational number null predicate. [q] is a RAT. Returns [t] or
% [nil]. Double null is allowed for.
null numr q or numr q = 0;
procedure rat_sgn(q);
% Rational number sign. [q] is of type RAT. Returns -1, 0 or 1
if null numr q or numr q = 0 then
0
else if sgn numr q = sgn denr q then
1
else
-1;
procedure rat_0;
% Rational number 0. Returns a RAT.
simp 0;
procedure rat_1;
% Rational number 1. Returns a RAT.
simp 1;
procedure rat_less(p,q);
% Rational number less. [p],[q] are RAT's. Returns [t] or [nil].
rat_sgn(addsq(p,negsq q)) = -1;
procedure rat_greater(p,q);
% Rational number greater. [p],[q] are RAT's. Returns [t] or [nil].
rat_sgn(addsq(p,negsq q)) = 1;
procedure rat_eq(p,q);
% Rational number equal. [p], [q] are RAT's. Returns [t] or [nil].
rat_sgn(rat_minus(p,q)) eq 0;
procedure rat_leq(p,q);
% Rational number less equal. [p], [q] are RAT's. Returns [t] or
% [nil].
rat_sgn(addsq(p,negsq q)) < 1;
procedure rat_min(p,q);
% Rational number minimum. [p], [q] are RAT's. Returns a RAT.
if rat_leq(p,q) then
p
else
q;
procedure rat_max(p,q);
% Rational number minimum. [p], [q] are RAT's. Returns a RAT.
if rat_leq(p,q) then
q
else
p;
procedure rat_mapmax(pl);
% Rational number map maximum. [pl] is a non-empty list of RAT.
% Returns a RAT.
if null cdr pl then
car pl
else
rat_max(car pl,rat_mapmax cdr pl);
procedure rat_min4(a,b,c,d);
% Rational number minimum of 4. [a], [b], [c], [d] are
% RAT's. Returns a RAT.
rat_min(rat_min(a,b),rat_min(c,d));
procedure rat_max4(a,b,c,d);
% Rational number maximum of 4. [a], [b], [c], [d] are
% RAT's. Returns a RAT.
rat_max(rat_max(a,b),rat_max(c,d));
procedure rat_abs(p);
% Rational number absolut value. [p] is a RAT. Returns a RAT.
absf numr p ./ denr p;
% module ratpoly;
%DS
% <RATPOLY> ::= <SQ> over integer, where the denominator is an integer
procedure ratpoly_fromatom(a);
% Rational number polynomial atom to rational number
% polynomial. [a] is an atom. Returns a RATPOLY.
simp a;
procedure ratpoly_fromsf(f);
% Rational number polynomial standard form to rational number
% polynomial. [f] is a SF over integers. Returns an RATPOLY.
!*f2q f;
procedure ratpoly_fromrat(r);
% Rational number polynomial from rational number. [r] is a
% RAT. Returns an RATPOLY.
r;
procedure ratpoly_torat(r);
<< if !*rlanuexdebug and ratpoly_idl r then
prin2t "***** ratpoly_torat: argument is not a rational";
r
>>;
procedure ratpoly_toaf(r);
prepsq r;
procedure ratpoly_0();
simp 0;
procedure ratpoly_1();
simp 1; % ratpoly_fromatom(1);
procedure ratpoly_mklin(x,ba);
% makes linear polynomial a*x-b, which has the root ba.
ratpoly_fromsf addf(multf(numr simp x,rat_denr ba),negf rat_numr ba);
procedure ratpoly_print q;
mathprint prepsq q;
procedure ratpoly_ratp(q);
% Rational number polynomial rational number predicate. [q] is a
% RATPOLY. Returns [t] or [nil]. Tests, if a RATPOLY is a RAT.
numberp numr q or null numr q;
procedure ratpoly_null();
% Rational number polynomial null. Returns a RATPOLY.
simp 0;
procedure ratpoly_nullp(q);
% Rational number polynomial null predicate. [q] is a
% RATPOLY. Returns [t] or [nil].
null numr q or numr q = 0;
procedure ratpoly_sgn q;
%%% ratpolydebug
if !*rlanuexdebug and not ratpoly_ratp q then
prin2t "***** ratpoly_sgn: not a constant polynomial"
else
rat_sgn q;
procedure ratpoly_neg(q);
% Rational number polynomial negation. [q] is a RATPOLY. Returns a
% RATPOLY.
negsq(q);
procedure ratpoly_add(q1,q2);
% Rational number polynomial addition. [q1], [q2] are
% RATPOLY's. Returns a RATPOLY.
addsq(q1,q2);
procedure ratpoly_minus(q1,q2);
% Rational number polynomial minus. [q1], [q2] are
% RATPOLY's. Returns a RATPOLY.
subtrsq(q1,q2);
procedure ratpoly_mult(q1,q2);
% Rational number polynomial multiplication. [q1], [q2] are
% RATPOLY's. Returns a RATPOLY.
multsq(q1,q2);
procedure ratpoly_foldmult(ql);
if null ql then
ratpoly_1()
else
ratpoly_mult(car ql,ratpoly_foldmult(cdr ql));
procedure ratpoly_quot(q1,q2);
% Rational number polynomial quotient. [q1], [q2] are
% RATPOLY's, [q2] is not null. Returns a RATPOLY.
quotsq(q1,q2);
procedure ratpoly_pp(q);
(sfto_dprpartksf numr q) . denr q;
procedure sf_univarp(f);
if domainp f or (domainp lc f and sf_univarp1(red f,mvar f)) then
t;
procedure sf_univarp1(f,x);
if domainp f or (domainp lc f and mvar f eq x and sf_univarp1(red f,x)) then
t;
procedure ratpoly_univarp(q);
sf_univarp numr q;
procedure ratpoly_idl q;
% Identifier list.
sf_idl numr q;
procedure ratpoly_mvartest(q,x);
% Rational number polynomial main variable test. [q] is a RATPOLY,
% [x] is an identifier. Returns [t] or [nil].
if null numr q then
nil
else if numberp numr q then
nil
else
(mvar numr q eq x);
procedure ratpoly_lc(q);
% Rational number polynomial leading coefficient. [q] is a RATPOLY.
% Returns a RATPOLY, the leading coefficient of $q$.
numpoly_lc numr q ./ denr q;
procedure ratpoly_red(q);
% Rational number polynomial reductum. [q] is a a RATPOLY.
% Returns a RATPOLY, the reductum of $q$.
numpoly_red numr q ./ denr q;
procedure ratpoly_ldeg(q);
% Rational number polynomial leading degree. [q] is a
% RATPOLY. Returns -1 or a natural number.
numpoly_ldeg numr q;
procedure ratpoly_deg(q,x);
% Degree of [q] in [x].
if ratpoly_mvartest(q,x) then
ratpoly_ldeg(q)
else
if ratpoly_nullp q then -1 else 0;
procedure ratpoly_exp(rp,n);
if n eq 0 then
ratpoly_1()
else
ratpoly_mult(rp,ratpoly_exp(rp,n-1));
procedure ratpoly_xtothen(x,n);
% Rational number polynomial x to the n. [x] is an identifier, [n]
% is a natural number. Returns a RATPOLY.
if n = 0 then
ratpoly_fromatom 1
else
ratpoly_mult(ratpoly_fromatom x,ratpoly_xtothen(x,n-1));
procedure ratpoly_diff(q,x);
% Rational number polynomial derivation. [q] is a RATPOLY, [x] is
% an identifier. Returns a RATPOLY.
diffsq(q,x);
procedure ratpoly_subrat(q,k,r);
% Rational number polynomial substitute rational number. [q] is a
% RATPOLY, [k] an identifier (the kernel to be replaced) and [r] is
% a SQ (the replacement). Returns a RATPOLY, the exact result of
% the substitution.
if domainp numr q or mvar numr q neq k then q else
% this might be faster here. test, if needed.
% !*f2q(cdr qremf(multf(numr q,exptf(denr r,ldeg numr q)),
% addf(multf(!*k2f k,denr r),negf numr r)));
begin scalar cl,res;
cl := coeffs numr q; % (c0,...,cn)
res := !*f2q car cl;
for each c in cdr cl do
res := addsq(!*f2q c,multsq(res,r));
res := quotsq(res,!*f2q denr q);
if !*rlanuexdebug and not ratpoly_nullp
ratpoly_minus(res,subsq(q,{k . prepsq r})) then
prin2 "***** ratpoly_subrat: faulty calculation";
return res;
end;
procedure ratpoly_subrat1(q,k,r);
% Rational number polynomial substitute rational number. [q] is a
% RATPOLY, [k] an identifier (the kernel to be replaced) and [r] is
% a SQ (the replacement). Returns a RATPOLY. Remark: if opt is nil
% then the the result is exact, otherwise, if opt is t, at least
% the sign is correct. %subsq(q,{k . prepsq r});
if domainp numr q or mvar numr q neq k then q else
begin scalar cl,res;
cl := coeffs numr q; % (c0,...,cn)
res := !*f2q car cl;
for each c in cdr cl do
res := addsq(!*f2q c,multsq(res,r));
if !*rlanuexdebug and not ratpoly_nullp
ratpoly_minus(quotsq(res,!*f2q denr q),
subsq(q,{k . prepsq r})) then
prin2 "***** ratpoly_subrat: faulty calculation";
return res;
end;
procedure ratpoly_sub(q,k,af);
% Rational number polynomial substitute rational number. [q] is a
% RATPOLY, [k] an identifier (the kernel to be replaced) and [af]
% is an algebraic form (the replacement). Returns a RATPOLY.
subsq(q,{k . af});
procedure ratpoly_tad(q);
% Throw away denominator.
numr q ./ 1;
procedure ratpoly_gcd(q1,q2); % plus
% Rational number polynomial gretest common divisor. [q1], [q2] are
% RATPOLY's. Returns a RATPOLY.
begin scalar tmp1,tmp2,res;
% convert to algebraic expression
tmp1 := prepsq q1;
tmp2 := prepsq q2;
on1 'rational;
% convert to sf over rationals
tmp1 := numr simp tmp1;
tmp2 := numr simp tmp2;
res := gcdf(tmp1,tmp2);
% convert to algebraic expression
res := prepf res;
off1 'rational;
return simp res;
end;
procedure ratpoly_resultant(q1,q2,y);
% Rational number polynomial gretest common divisor. [q1], [q2] are
% RATPOLY's. Returns a RATPOLY.
begin scalar tmp1,tmp2,res;
% convert to algebraic expression
tmp1 := prepsq q1;
tmp2 := prepsq q2;
on1 'rational;
% convert to sf over rationals
tmp1 := numr simp tmp1;
tmp2 := numr simp tmp2;
res := resultant(tmp1,tmp2,y);
% convert to algebraic expression
res := prepf res;
off1 'rational;
return simp res;
end;
procedure ratpoly_factorize(q);
% Type: ratpoly -> list(pair(ratpoly,num))
begin scalar tmp;
tmp := numpoly_factorize numr q; % gives a
% const factor has multiplicity 1
car tmp := ratpoly_fromrat rat_mk(car tmp,denr q) . 1;
tmp := car tmp .
for each sfn in cdr tmp collect (ratpoly_fromsf car sfn) . cdr sfn;
if !*rlanuexdebug and not
ratpoly_nullp(ratpoly_minus(q,ratpoly_foldmult(
for each sfn in tmp collect ratpoly_exp(car sfn,cdr sfn)))) then
prin2t "***** ratpoly_factorize: selftest failed";
return tmp;
end;
% module sf;
procedure sf_deg(f,x);
if domainp f then
if null f then -1 else 0
else
if mvar f = x then ldeg f else 0;
% module ctx;
% Context
% IAL ::= LIST(PAIR(ID,ANU))
% Type: {'ctx, ial} where ial
% Description:
procedure ctx_fromial(ia);
{'ctx,ia};
procedure ctx_new;
ctx_fromial nil;
procedure ctx_ial(c);
cadr c;
procedure ctx_idl(c);
for each ia in ctx_ial c collect car ia;
procedure ctx_print(c);
<<
prin2 "[(ctx) ";
for each ia in ctx_ial c do
<<prin2 car ia; prin2 "->"; anu_print cdr ia>>;
prin2 "] "
>>;
procedure ctx_add(ia,c);
% Add, no check. [ia] is a dotted pair (xj . a) where xj is an id
% an a is an anu.
{'ctx,ia . ctx_ial c};
procedure ctx_addip(ia,c);
% Add, no check. [ia] is a dotted pair (xj . a) where xj is an id
% an a is an anu.
cadr c := ia . ctx_ial c;
procedure ctx_rm(x,c);
ctx_fromial(for each ia in ctx_ial c join
if car ia eq x then nil else {ia});
procedure ctx_free(x,c);
% frees x and all higher variables. x has to be bound by c.
begin scalar idorder,ial;
if !*rlanuexdebug and null ctx_get(x,c) then
prin2t "***** ctx_free: variable unbound";
idorder := setkorder nil;
setkorder idorder;
ial := ctx_ial c;
repeat <<
while caar ial neq car idorder do idorder := cdr idorder;
ial := cdr ial; % free the highest pair
>> until null idorder or car idorder eq x;
return ctx_fromial ial;
end;
procedure ctx_get(x,c);
% Returns the algebraic number to which x is bound, nil if x is
% unbound.
begin scalar ial,res;
if !*rlanuexdebug and not idp x then rederr "ctx_get: arguments invalid";
ial := ctx_ial c;
while ial and null res do <<
if caar ial eq x then res := cdar ial; ial := cdr ial
>>;
if !*rlanuexdebug and null res then
prin2 "***** ctx_get: variable unbound";
return res;
end;
procedure ctx_union(c1,c2);
% Union of syntactically compatible contexts.
% Syntactically/semantically compatible means: for every id x: if
% (x,a1) in c1 and (x,a2) in c2 then syntactically/semantically
% a1=a2.
begin scalar idorder,ial;
idorder := setkorder nil;
setkorder idorder;
c1 := ctx_ial c1; c2 := ctx_ial c2;
for each id in idorder do %%% while loop might be faster
if not null c1 and caar c1 eq id then <<
ial := car c1 . ial;
c1 := cdr c1;
if c2 then if caar c2 eq id then c2 := cdr c2;
>>
else if not null c2 and caar c2 eq id then <<
ial := car c2 . ial;
c2 := cdr c2;
>>;
if c1 or c2 then prin2t "***** ctx_union: idorder not complete";
return ctx_fromial reversip ial;
end;
% module aex;
% Algebraic expressions.
%%% --- constructors and access functions --- %%%
procedure aex_fromrp rp;
{'aex,rp,ctx_new(),t,t};
procedure aex_fromrat r;
{'aex,ratpoly_fromrat r,ctx_new(),t,t};
procedure aex_0();
aex_fromrp ratpoly_0();
procedure aex_1();
aex_fromrp ratpoly_1();
procedure aex_mk(rp,c,lcnttag,reducedtag);
% Make. [rp] is a ratpoly, [c] is a context. Return an aex.
{'aex,rp,c,lcnttag,reducedtag};
procedure aex_mklin(x,ba);
% Make linear poly a*x+b. x is an id, ba is a rat.
%aex_fromsf addf(multf(numr simp x,rat_denr ba),negf rat_numr ba);
aex_fromrp ratpoly_mklin(x,ba);
procedure aex_fromsf(f);
aex_fromrp ratpoly_fromsf f;
procedure aex_mapfromsf(fl);
for each f in fl collect aex_fromsf f;
procedure aex_ex(ae);
cadr ae;
procedure aex_ctx(ae);
caddr ae;
procedure aex_freeall(ae);
aex_mk(aex_ex ae,ctx_new(),t,t);
procedure aex_free(ae,x);
% Free x and all higher variables.
aex_mk(aex_ex ae,ctx_free(x,aex_ctx ae),nil,nil);
procedure aex_bind(ae,x,a);
%%% todo: ensure [a] is defined with x
% test if [a] is rational (then use aex_subrat)
if null aex_boundids anu_dp a and aex_deg(anu_dp a,x) eq 1 then
aex_subrp(ae,x,ratpoly_toaf aex_ex aex_linsolv(anu_dp a,x))
else aex_mk(aex_ex ae,ctx_add(x . a,aex_ctx ae),nil,nil);
procedure aex_print(ae);
<<
ratpoly_print aex_ex ae;
if ctx_ial aex_ctx ae then << prin2 ", where ";ctx_print aex_ctx ae; >>
>>;
% - tags stuff - %
procedure aex_lcnttag ae;
cadddr ae;
procedure aex_reducedtag ae;
cadr cdddr ae;
procedure aex_putlcnttag(ae,v);
cadddr ae := v;
procedure aex_putreducedtag(ae,v);
cadr cdddr ae := v;
%% --- arithmetic with pseudodivision, sign, nullp, psgcd pssqfree --- %%%
procedure aex_neg(ae);
aex_mk(ratpoly_neg aex_ex ae,aex_ctx ae,aex_lcnttag ae,aex_reducedtag ae);
procedure aex_add(ae1,ae2);
% Add, contexts are assumed to be compatible and will be merged.
% caveat: minimization will be needed.
aex_mk(ratpoly_add(aex_ex ae1,aex_ex ae2),
ctx_union(aex_ctx ae1, aex_ctx ae2),
nil,aex_reducedtag ae1 and aex_reducedtag ae2);
procedure aex_addrat(ae,r);
% Add rational number.
aex_mk(ratpoly_add(aex_ex ae, ratpoly_fromrat r),aex_ctx ae,
nil,nil); %%%?
procedure aex_minus(ae1,ae2);
% Minus, contexts are assumed to be compatible and will be merged.
aex_mk(ratpoly_minus(aex_ex ae1,aex_ex ae2),
ctx_union(aex_ctx ae1, aex_ctx ae2),
nil,aex_reducedtag ae1 and aex_reducedtag ae2);
procedure aex_mult(ae1,ae2);
% Multiplication, contexts are assumed to be compatible and will be merged.
aex_mk(ratpoly_mult(aex_ex ae1,aex_ex ae2),
ctx_union(aex_ctx ae1, aex_ctx ae2),
aex_lcnttag ae1 and aex_lcnttag ae2,nil);
procedure aex_foldmult(ael);
if null ael then
aex_1()
else
aex_mult(car ael,aex_foldmult cdr ael);
procedure aex_multrat(ae,r);
% Multiply with rational number.
aex_mk(ratpoly_mult(aex_ex ae, ratpoly_fromrat r),aex_ctx ae,
nil,nil);
procedure aex_diff(ae,x);
% differentiate. [ae] is an algebraic polynomial in [x].
% Returns an AEX.
aex_mk(ratpoly_diff(aex_ex ae,x),aex_ctx ae,
nil,nil);
procedure aex_subrp(ae,x,af);
% Substitute algebraic form in algebraic expression.
aex_mk(ratpoly_sub(aex_ex ae,x,af),aex_ctx ae,
nil,nil); %%% minimize ex and ctx
procedure aex_subrat(ae,x,r);
% exact
aex_mk(ratpoly_subrat(aex_ex ae,x,r),aex_ctx ae,nil,nil);
procedure aex_subrat1(ae,x,r);
% exact up to sign
aex_mk(ratpoly_subrat1(aex_ex ae,x,r),aex_ctx ae,nil,nil);
procedure aex_tad(ae);
% Throw away denominator.
aex_mk(ratpoly_tad aex_ex ae,aex_ctx ae,nil,nil);
procedure aex_xtothen(x,n);
% x to the n.
aex_mk(ratpoly_xtothen(x,n),ctx_new(),
t,t);
procedure aex_deg(ae,x);
% Degree of [ae] in [x].
ratpoly_deg(aex_ex ae,x);
procedure aex_simpleratpolyp(ae);
% Simple but uncomplete predicte to test, if ae represents a
% (maybe constant) rational polynomial
null ctx_ial aex_ctx ae or ratpoly_ratp aex_ex ae;
procedure aex_simpleratp ae;
ratpoly_ratp aex_ex ae;
procedure aex_simplenullp(ae);
% Simple but uncomplete null-predicate. It is complete, if ae is
% minimized.
ratpoly_nullp(aex_ex ae);
procedure aex_simplenumberp(ae);
not aex_freeids ae;
procedure aex_ids(ae);
ratpoly_idl aex_ex ae;
procedure aex_freeids(ae);
% free identifiers, the id with highest kernel order first
setminus(ratpoly_idl aex_ex ae,ctx_idl aex_ctx ae);
procedure aex_boundids(ae);
intersection(aex_ids ae,ctx_idl aex_ctx ae);
procedure aex_constp(ae);
% Constant predicate. %%% faster!
null aex_freeids ae;
procedure aex_nullp(ae);
% Null predicate. [ae] is an AEX. Returns [t] or [nil].
begin scalar tmp;
% make the lc non-trivial
tmp := aex_mklcnt ae;
% ae is a non-constant polynomial
if aex_freeids tmp then
return nil;
% ae is a constant
if aex_sgn tmp eq 0 then
return t
else
return nil;
end;
procedure aex_mvaroccurtest(ae,x);
ratpoly_mvartest(aex_ex ae,x);
procedure aex_red(ae,x);
% Reductum of [ae] wrt [x]. Needs not to be minimized.
if aex_mvaroccurtest(ae,x) then
aex_mk(ratpoly_red aex_ex ae,aex_ctx ae,nil,nil) %%% mklcnt
else
aex_0();
procedure aex_lc(ae,x);
if aex_mvaroccurtest(ae,x) then
aex_mk(ratpoly_lc aex_ex ae,aex_ctx ae,nil,nil)
else
ae; % ctx needs not to be made smaller, as there are no singles
procedure aex_mvar ae;
begin scalar idl;
% semanticcheck
if not (idl := aex_freeids ae) then
prin2t "***** aex_mvar: invalid argument";
return car idl;
end;
procedure aex_mklcnt(ae);
% Make leading coefficient of non-constant polynomial non-trivial.
% was: minimize.
begin scalar idl;
% quick win: obvious cases: ae is rat or contains no alg. numbers
%% turn around, use ctx_ial
if ratpoly_ratp aex_ex ae or null ctx_idl aex_ctx ae then
return ae;
% ae is a non-constant algebraic polynomial
if idl := aex_freeids ae then
if aex_nullp aex_lc(ae,car idl) then
return aex_mklcnt aex_red(ae,car idl)
else
return ae;
% ae is constant
return %%% remark: this is the second edge in the dependence graph where the argument is passed and not made smaller. (this is not a problem)
if aex_sgn ae eq 0 then aex_0() else ae
end;
%procedure aex_minimize(ae);
% % Algebraic expression minimize. [ae] is an AEX(c,d) with
% % $c<d$. Returns an AEX(c,d), where the $r$-component is null or
% % where the leading coefficient is non-trivial.
% if aex_getr ae = ratpoly_null() then
% aex_mkfromr(ae,ratpoly_null())
% else if aex_nullp aex_lc ae then
% aex_minimize aex_red ae
% else
% ae;
%procedure aex_reduce(ae);
% % convention: the bound variable has to be defined by an anu using
% % the same variable.
% begin scalar ids,x,alpha,rlc,rred,tmp;
% %%% care for special case aex_simpleratp !!!
% % there are no bound variables: we have a (const.) rat. poly
% if not aex_boundids ae then
% return ae;
% % from now on there are bound variables
% ids := aex_ids ae; x := car ids;
% % there are free variables
% if aex_freeids ae then <<
% rlc := aex_reduce(aex_lc(ae,x));
% rred := aex_reduce(aex_red(ae,x));
% return aex_add(aex_mult(rlc,aex_xtothen(x,aex_deg(ae,x))),rred)
% >>;
% % there are no free variables
% if !*rlanuexdebug and not (x eq caar ctx_ial aex_ctx ae) then
% prin2t "aex_reduce: variable mismatch";
% %%%alpha := cdar ctx_ial aex_ctx ae; % context = {(x . alpha),...}
% alpha := ctx_get(x,aex_ctx ae);
% tmp := aex_free(ae,x);
% tmp := aex_reduce tmp;
% if aex_deg(tmp,x) >= aex_deg(anu_dp alpha,x) then
% return aex_bind(aex_rem(tmp,anu_dp alpha,x),x,alpha)
% else
% return aex_bind(tmp,x,alpha)
% end;
procedure aex_reduce ae;
% convention: the bound variable has to be defined by an anu using
% the same variable.
begin scalar ids,x,alpha,rlc,rred,tmp;
% there are no bound variables: we have a (const.) rat. poly
if not aex_boundids ae then
tmp := ae;
% from now on there are bound variables
if null tmp then << ids := aex_ids ae; x := car ids >>;
% there are free variables
if null tmp and aex_freeids ae then <<
rlc := aex_reduce(aex_lc(ae,x));
rred := aex_reduce(aex_red(ae,x));
tmp := aex_add(aex_mult(rlc,aex_xtothen(x,aex_deg(ae,x))),rred)
>>
% there are no free variables
else if null tmp then <<
%if !*rlanuexdebug and not (x eq caar ctx_ial aex_ctx ae) then
% prin2t "aex_reduce: variable mismatch";
%alpha := cdar ctx_ial aex_ctx ae; % context = {(x . alpha),...}
alpha := ctx_get(x,aex_ctx ae);
tmp := aex_free(ae,x);
tmp := aex_reduce tmp;
if aex_deg(tmp,x) >= aex_deg(anu_dp alpha,x) then
tmp := aex_bind(aex_rem(tmp,anu_dp alpha,x),x,alpha)
else
tmp := aex_bind(tmp,x,alpha)
>>;
% if !*rlanuexdebug and not aex_nullp aex_minus(ae,tmp) then
% prin2t "aex_reduce: selftest failed";
return tmp;
end;
procedure aex_reducetest(ae);
aex_nullp aex_minus(ae,aex_reduce ae);
procedure aex_psquotrem1(f,p,x);
% pseudo quotient remainder one step. [f], [p] are algebraic
% polynomials in [x]. Returns [(q . r)] with $a_n^2f=qp+r$.
begin scalar am,an,m,n,q;
am := aex_lc(f,x); % a constant algebraic poly
an := aex_lc(p,x); % dito
m := aex_deg (f,x);
n := aex_deg (p,x);
if m<n or n=-1 then
return (aex_0(). f);
% 0<=m<=n
q := aex_mult(aex_xtothen(x,m-n),aex_mult(am,an));
return (q . aex_minus(aex_mult(f,aex_mult(an,an)),aex_mult(p,q)));
% return (q . aex_minus(aex_mult(aex_red f,aex_mult(an,an)),
% aex_red aex_mult(p,q)));
end;
procedure aex_psquotrem(f,p,x);
% pseudo quotient remainder. [f], [p] are algebraic polynomials in
% [x]. Returns [(q . r . i)] with $b_n^{2*i}f=qp+r$.
begin scalar m,n,q,r,bn2,qr1; integer i,ii;
m := aex_deg(f,x);
n := aex_deg(p,x);
q := aex_0();
r := f;
if m<n or aex_simplenullp(p) then
return (q . r . 0);
bn2 := aex_mult(bn,bn) % bn2 needed below for selftest
%bn2 := aex_multrat(bn,rat_fromnum aex_sgn bn) % wrong; change in psquotrem1 needed as well!
where bn = aex_lc(p,x);
while (aex_deg(r,x) >= n) do <<
qr1 := aex_psquotrem1(r,p,x);
q := aex_add(aex_mult(q,bn2),car qr1);
i := i+1;
r := cdr qr1;
>>;
if !*rlanuexdebug then <<
ii := i;
while (ii:=ii-1) >= 0 do % ii times multply bn^2
f := aex_mult(f,bn2);
if not aex_nullp(aex_minus(f,aex_add(aex_mult(q,p),r))) then
prin2t "***** aex_psquotrem: selftest failed";
>>;
return (q . r . i);
end;
procedure aex_psquot(f,p,x);
car aex_psquotrem(f,p,x);
procedure aex_psrem(f,p,x);
cadr aex_psquotrem(f,p,x);
procedure aex_psquotremtest(f,p,x);
% Algebraic expression pseudo quotient remainder test.
% Returns a BOOL.
begin scalar m,n,ctx,q,r,bn2,qri,ii;
m := aex_deg(f,x); n := aex_deg(p,x);
ctx := aex_ctx f;
qri := aex_psquotrem(f,p,x);
q := car qri;
r := cadr qri;
ii := cddr qri;
bn2 := aex_mult(bn,bn) where bn=aex_lc(p,x);
while (ii:=ii-1) >= 0 do % ii times multply bn^2
f := aex_mult(f,bn2);
if not aex_nullp(aex_minus(f,aex_add(aex_mult(q,p),r))) then
prin2 "***** aex_psquotremtest: failed";
end;
procedure aex_stdsturmchain(f,x);
% standard sturm chain. [f] is an algebraic polynomial in [x].
% Returns a list of algebraic polynomials.
aex_sturmchain(f,aex_diff(f,x),x);
procedure aex_sturmchain(f,g,x);
% pseudo sturm chain. [f], [g] are algebraic expressions in [x].
% Returns a list of algebraic expressions.
begin scalar sc;
sc := reversip(aex_remseq({ aex_tad g,aex_tad f},x));
%%% sc := reversip(aex_remseq({aex_pp aex_tad g,aex_pp aex_tad f},x));
%if !*rlanuexdebug then aex_sturmchaincheck sc;
return sc
end;
procedure aex_sturmchaincheck(sc);
% Algebraic expression sturm chain check. [sc] is a sturm
% chain. Returns a BOOL, and prints an error message.
begin scalar v,l;
v := t;
l := length sc;
while sc and v do
if aex_nullp car sc then <<
v := nil;
ioto_prin2 {"+++ aex_sturmchaincheck: nullpoly found. length:",
l," position:",l-length sc+1," +++"}
>> else
sc := cdr sc;
return v
end;
procedure aex_pp(ae);
% primitive part. works only for univariate polynomials with
% rational coefficients.
% ae;
<< if !*rlanuexdebug and aex_simpleratpolyp ae and
length aex_ids ae > 1 then
prin2 "***** aex_pp: argument not univariate";
if ratpoly_univarp aex_ex ae then
aex_mk(ratpoly_pp aex_ex ae,aex_ctx ae,nil,nil)
else
ae
>>;
procedure aex_remseq(ael,x);
% remainder sequence. [ael] is a list of algebraic polynomials in
% [x], of at least length 2. Returns a list of algebraic
% polynomials. Caveat: the list is built in reverse order.
begin scalar rem;
if !*rlanuexdebug and aex_simplenullp car ael then
prin2 "[remseq:null]";
if aex_deg(car ael,x) <= 0 then
return ael;
if !*rlanuexpsremseq then
rem := aex_psrem(cadr ael,car ael,x)
else
rem := aex_pp aex_tad aex_rem(cadr ael,car ael,x);
if aex_simplenullp rem then
return ael;
return aex_remseq(aex_neg rem . ael,x)
end;
procedure aex_sturmchainsgnch(sc,x,r);
% Algebraic expression sturm chain sign changes. [sc] is a list of
% AEX(c,c+1), [r] is a RAT. Returns n integer, the number of sign
% changes of $sc$ an $r$.
num_sgnch for each e in sc collect aex_sgn aex_subrat1(e,x,r);
procedure aex_stchsgnch(sc,x,r);
% sturm chain sign changes extended version. [r] is a rat or infty
% or minfty.
if r eq 'infty then
num_sgnch for each e in sc collect aex_sgnatinfty(e,x)
else if r eq 'minfty then
num_sgnch for each e in sc collect aex_sgnatminfty(e,x)
else
aex_sturmchainsgnch(sc,x,r);
procedure aex_sgnatinfty(ae,x);
% Sign at infinity. [ae] has non-trivial lc or is simply null.
begin scalar freeids;
if aex_simplenullp ae then return 0;
freeids := aex_freeids ae;
if not freeids then return aex_sgn ae;
% from now on we have one free variable.
if !*rlanuexdebug and car freeids neq x then
prin2 "***** aex_sgnatinfty: wrong main variable";
if !*rlanuexdebug and length freeids > 1 then
prin2 "***** aex_sgnatinfty: not univariate";
return aex_sgn aex_lc(ae,x);
end;
procedure aex_sgnatminfty(ae,x);
% Sign at minus infinity. [ae] has non-trivial lc or is simply null.
begin scalar freeids;
if aex_simplenullp ae then return 0;
freeids := aex_freeids ae;
if not freeids then return aex_sgn ae;
% from now on we have one free variable
if !*rlanuexdebug and car freeids neq x then
prin2 "***** aex_sgnatminfty: wrong main variable";
if !*rlanuexdebug and length freeids > 1 then
prin2 "***** aex_sgnatminfty: not univariate";
if num_even aex_deg(ae,x) then
return aex_sgn aex_lc(ae,x);
return (-1)*aex_sgn aex_lc(ae,x);
end;
procedure aex_sgn(ae);
% Remark: the case: ae has free ids, but is constant due to trivial
% coefficients is not treated, although sgn would be sensible.
% Possible optimization: use of aex_containment.
begin scalar con,x,g,alpha,f,sc;
% semanticcheck
if !*rlanuexdebug and (not (type_of ae eq 'aex) or aex_freeids ae) then
prin2t "***** aex_sgn: invalid argument";
% ae is obviously rational %%% faster with _ratp
if ratpoly_ratp aex_ex ae then
return ratpoly_sgn aex_ex ae;
% possible optimization:
if !*rlanuexsgnopt then <<
con := aex_containment ae;
if rat_less(rat_0(),iv_lb con) then return 1;
if rat_less(iv_rb con,rat_0()) then return (-1);
>>;
% ae is algebraic and ae=(g(x),(x=f(y),etc))
%x := caar ial; % flawed, if ctx is not minimal
x := car ratpoly_idl aex_ex ae;
% alpha := cdar ial; % flawed, if ctx is not minimal
alpha := ctx_get(x,aex_ctx ae);
g := aex_mklcnt aex_reduce aex_free(ae,x); %%% reduce somewhere else, eg bind
if aex_simpleratp g then return ratpoly_sgn aex_ex g;
if ofsf_anuexverbosep() and aex_deg(g,x)<=0 then prin2 "[aex_sgn:num!]";
f := anu_dp alpha;
f := aex_subrp(f,aex_mvar f,x); %unnecessary, aex_bind makes it already
if !*rlanuexdebug and aex_sgn aex_lc(f,x) eq 0 then
prin2t "***** aex_sgn: f has trivial lc";
if !*rlanuexdebug and aex_sgn aex_lc(g,x) eq 0 then
prin2t "***** aex_sgn: g has trivial lc";
% sc is the sturmchain for f(x), f'g(x)
%%% optimization: aex_reduce after aex_diff
sc := aex_sturmchain(f,aex_mult(aex_diff(f,x),g),x);
return aex_sturmchainsgnch(sc,x,iv_lb anu_iv alpha)
- aex_sturmchainsgnch(sc,x,iv_rb anu_iv alpha)
end;
procedure aex_pssqfree(f,x);
% pseudo sqare-free part. [f] is an algebraic polynomial. Returns
% an algebraic polynomial.
if aex_deg(f,x) < 2 then f
else car aex_psquotrem(f,lastcar aex_stdsturmchain(f,x),x);
%procedure aex_psgcd(f,g,x);
% % pseudo greatest common divisor. [f], [g] are algebraic
% % polynomials with positive degree in [x]. Returns an algebraic
% % polynomial.
% %lastcar aex_sturmchain(f,g);
% car aex_remseq({g,f},x);
%%% --- should be after exact arithmetic --- %%%
procedure sqfr_norm(f,x,y,palpha);
% (RATPOLY,ID,ID,RATPOLY)->(NUM,RATPOLY,RATPOLY)
begin scalar g,r; integer s,degree;
s := 0; g := f;
palpha := ratpoly_sub(palpha,x,y);
repeat <<
r := ratpoly_resultant(palpha,g,y);
% ratpoly_deg works just if x is the leading kernel
degree := ratpoly_deg(ratpoly_gcd(r,ratpoly_diff(r,x)),x);
if degree neq 0 then <<
s := s+1;
g := ratpoly_sub(g,x,{'plus,x,{'minus,y}})
>>
>>
until degree = 0;
return {s,g,r};
end;
procedure aex_factorize(f,x); %%% rename to factors
% Factors. [f] is a squarefree univariate alg. polynomial in [x].
% Returns the list of non-constant factors of [f]. Remark: this
% implements trager's alg_factor. Funthermore, if
% [aex_simpleratpolyp] holds, [f] needs not to be squarefree.
begin scalar tmp,y,alpha,s,g,r,l,aeg,aehi;
if aex_simpleratpolyp f then <<
l := ratpoly_factorize aex_ex f;
% self-test done already in ratpoly_factorize.
return cdr for each rp_m in l collect aex_fromrp car rp_m;
>>;
tmp := car ctx_ial aex_ctx f; % y . alpha
y := car tmp; alpha := cdr tmp;
%tmp := aex_sqfr_norm(f,x,y); % maybe in future...
tmp := sqfr_norm(aex_ex f,x,y,aex_ex anu_dp alpha);
s := car tmp; g := cadr tmp; r := caddr tmp;
tmp := ratpoly_factorize r; %%%%%%%
% throw away multiplicities and domain element
l := cdr for each rp_m in tmp collect car rp_m;
%l := for each rp_m in l collect car rp_m;
aeg := aex_fromrp g; aeg := aex_bind(aeg,y,alpha); % g(x,alpha)
l := for each hi in l collect <<
aehi := aex_fromrp hi; % h_i(x)
aehi := aex_gcd(aehi,aeg,x); % h_i(x,alpha)=gcd(h_i(x),g(x,alpha))
aeg := aex_quot(aeg,aehi,x);
aehi := aex_subrp(aehi,x,{'plus,x,{'times,s,y}})
>>;
if !*rlanuexdebug and not
aex_simpleratp aex_quot(f,aex_foldmult l,x) then
prin2t "***** aex_factorize: selftest failed (alg. case)";
return l; %%% minimize every element first?
end;
%%% --- exact arithmetic - with inv --- %%%
%procedure aex_quotrem(f,g,x);
% % optimization: g is constant.
% if null aex_freeids g then aex_mult(aex_inv g,f) . aex_0()
% else aex_quotrem1(f,g,x);
procedure aex_quotrem(f,g,x);
% f,g need to have non-trivial leading coefficient. Literatur: see
% becker, weispfenning, kredel: groebner bases, p.81.
% property: rr is always returned with non-trivial lc.
begin scalar rr,gg,qq,an,bm,qqi; integer m,n;
if aex_simplenullp g then
prin2t "***** aex_quotrem: g=0"; %%% return here 0 . f
% optimization: g is constant.
if null aex_freeids g then
return aex_mult(aex_inv g,f) . aex_0();
if !*rlanuexdebug and aex_sgn aex_lc(f,x) eq 0 then
prin2t "***** aex_quotrem: first argument invalid";
if !*rlanuexdebug and aex_sgn aex_lc(g,x) eq 0 then
prin2t "***** aex_quotrem: second argument invalid";
rr := f; gg := g; qq := aex_0();
n := aex_deg(rr,x); m := aex_deg(gg,x);
while not aex_simplenullp rr and n >= m do <<
an := aex_lc(rr,x);
bm := aex_lc(gg,x);
qqi := aex_reduce aex_mult(aex_mult(an,aex_inv bm),
aex_xtothen(x,n-m)); % reduce added as optimization
rr := aex_mklcnt aex_minus(rr,aex_mult(qqi,gg)); %%% mit redukta
qq := aex_add(qq,qqi);
n := aex_deg(rr,x); m := aex_deg(gg,x);
>>;
if !*rlanuexdebug and not aex_nullp aex_minus(f,aex_add(aex_mult(qq,g),rr))
then prin2t "aex_quotrem: selftest failed";
return qq . rr
end;
procedure aex_quotremtest(f,g,x);
% should result in 0
aex_minus(f,aex_add(aex_mult(car qr,g),cdr qr))
where qr=aex_quotrem(f,g,x);
procedure aex_quot(f,g,x);
car aex_quotrem(f,g,x);
procedure aex_rem(f,g,x);
%if null aex_freeids g then aex_0() else cdr aex_quotrem1(f,g,x);
cdr aex_quotrem(f,g,x);
procedure aex_gcdext(a,b,x);
aex_gcdext1(a,b,x,!*rlanuexgcdnormalize);
procedure aex_gcdext1(a,b,x,gcdnormalize);
% extended euclidean algorithm, see becker, weispfenning, kredel:
% groebner bases, p.83.
begin scalar aa,bb,ss,tt,uu,vv,qr,ss1,tt1,tmp;
aa := a; bb := b;
ss := aex_1(); tt := aex_0();
uu := aex_0(); vv := aex_1();
while not aex_simplenullp bb do << %%% simplenullp should suffice
qr := aex_quotrem(aa,bb,x);
aa := bb; bb := cdr qr; % attention: cadr with psquotrem
ss1 := ss; tt1 := tt;
ss := uu; tt := vv;
uu := aex_minus(ss1,aex_mult(car qr,uu));
vv := aex_minus(tt1,aex_mult(car qr,vv));
>>;
if gcdnormalize then <<
tmp := aex_inv aex_lc(aa,x); aa := aex_mult(aa,tmp);
ss := aex_mult(ss,tmp); tt := aex_mult(tt,tmp);
>>;
if !*rlanuexdebug and not
aex_nullp aex_minus(aa,aex_add(aex_mult(ss,a),aex_mult(tt,b))) then
prin2t "***** aex_gcdext: selftest failed";
return {aa,ss,tt}
end;
procedure aex_gcd(a,b,x);
begin scalar d;
d := car aex_gcdext1(a,b,x,nil);
% optimizatfion: if the gcd is not a poly, then it's 1
if aex_deg(d,x)<1 then d := aex_1();
if !*rlanuexdebug and aex_nullp aex_lc(d,x) then
prin2 "*** aex_gcd: lc non-trivial";
return d;
end;
procedure aex_gcd1(a,b,x,gcdnormalize);
car aex_gcdext1(a,b,x,gcdnormalize);
procedure aex_gcdexttest(a,b,x);
% should result in 0
aex_minus(car ddsstt,
aex_add(aex_mult(cadr ddsstt,a),aex_mult(caddr ddsstt,b)))
where ddsstt = aex_gcdext(a,b,x);
procedure aex_pairwiseprime(ael,x);
% pairwise prime. ael is a list of AEX with lc non-trivial. Returns
% a list of AEX with lc non-trivial.
begin scalar pprestlist,tmp;
if length ael <= 1 then return ael;
pprestlist := aex_pairwiseprime(cdr ael,x);
tmp := aex_pairwiseprime1(car ael,pprestlist,x);
% lets assume aex_quot returns something with lc non-trivial.
return if not aex_simplenumberp tmp then
tmp . pprestlist else pprestlist
end;
procedure aex_pairwiseprime1(ae1,ae2l,x);
% makes ae1 prime to ae2l.
% << for each ae2 in ae2l do ae1 := aex_quot(ae1,aex_gcd(ae1,ae2,x),x);
% ae1
% >>;
<< while ae2l and not aex_simplenullp ae1 do <<
ae1 := aex_quot(ae1,aex_gcd(ae1,car ae2l,x),x);
ae2l := cdr ae2l;
>>; ae1 >>;
procedure aex_sqfree(f,x);
% pseudo sqare-free part. [f] is an algebraic polynomial. Returns
% an algebraic polynomial.
% question: what if f not reduced?
if aex_deg(f,x) < 2 then f
else car aex_quotrem(f,aex_gcd(f,aex_diff(f,x),x),x);
procedure aex_inv ae;
% invert a constant, not-null polynomial.
% ae represents an algebraic number, that is not null.
%%% the case of a purely rational ae occurs very often. avoid the progn by splitting into two functions.
begin scalar ia,x,alpha,g,f,d,f1,drs,tmp;
% ae is obviously a rational
if ratpoly_ratp aex_ex ae then
return aex_fromrp ratpoly_quot(ratpoly_1(),aex_ex ae);
% ae is an constant algebraic polynomial
ia := car ctx_ial aex_ctx ae; % (x . (anu f iv))
x := car ia;
alpha := cdr ia;
g := aex_free(ae,x);
f := anu_dp alpha;
d := car aex_gcdext(f,g,x);
%%% d := car aex_gcd(f,g,x); % why segmentation violation?
f1 := car aex_quotrem(f,d,x);
% apply ext. euclid. algorithm to f1,g, gives: 1=rf1+sg
drs := aex_gcdext(f1,g,x);
% caveat: car drs is non-zero rational, but i.g. <> 1
if !*rlanuexgcdnormalize then
tmp := aex_bind(aex_mult(car drs,caddr drs),x,alpha) %%% here alphanew with f1
else
tmp := aex_bind(aex_mult(aex_inv car drs,caddr drs),x,alpha);
if !*rlanuexdebug and
not aex_nullp aex_minus(aex_1(),aex_mult(ae,tmp)) then
rederr "aex_inv: selftest failed";
return tmp;
end;
procedure aex_invtest ae;
% should result in 0
aex_minus(aex_1(),aex_mult(ae,aex_inv ae));
procedure aex_linsolv(ae,x);
% [x] is the only free variable in ae.
<<
if !*rlanuexdebug and not (aex_freeids ae equal {x}) then
prin2t "aex_linsolv: invalid argument";
aex_mult(aex_inv aex_lc(ae,x),aex_neg aex_red(ae,x))
>>;
%%% --- root isolation --- %%%
procedure aex_coefl(ae,x);
if aex_mvaroccurtest(ae,x) then
aex_lc(ae,x) . aex_coefl(aex_red(ae,x),x)
else
{ae};
procedure aex_coefdegl(ae,x);
% coefficients and degree list. ae is a polynomial in x.
if aex_mvaroccurtest(ae,x) then
(aex_lc(ae,x) . aex_deg(ae,x)) . aex_coefdegl(aex_red(ae,x),x)
else
{aex_lc(ae,x) . 0};
procedure aex_fromcoefdegl(cfdgl,x);
begin scalar ae;
ae := aex_0();
for each cd in cfdgl do
ae := aex_add(ae,aex_mult(car cd,aex_xtothen(x,cdr cd)));
return ae
end;
procedure aex_coefdegltest(ae,x);
aex_nullp aex_minus(ae,aex_fromcoefdegl(aex_coefdegl(ae,x),x));
procedure aex_containment(ae);
% Algebraic expression containment. [ae] is an AEX. Returns an
% interval. the algebraic number represented by ae is contained in
% the interval, and the interval is regarded as closed.
begin scalar ia,cfdgl,ctac,ivl,r;
% coefficient and degree list, containment of a_c, interval list
if !*rlanuexdebug and aex_freeids ae then
prin2t "***** aex_containment: invalid argument";
if null aex_boundids ae then <<
r := ratpoly_torat aex_ex ae;
return iv_mk(r,r);
>>;
% there is a bound variable now
ia := car ctx_ial aex_ctx ae;
ctac := anu_iv cdr ia;
%ae := aex_free(ae,car ia);
cfdgl := aex_coefdegl(aex_free(ae,car ia),car ia);
if !*rlanuexdebug and not aex_coefdegltest(ae,car ia) then
prin2t "***** aex_containment: aex_coefdegltest failed";
ivl := for each cfdg in cfdgl collect
iv_mult(aex_containment car cfdg,iv_tothen(ctac,cdr cfdg));
return iv_mapadd ivl
end;
procedure aex_distinguishpositivefromzero(ae,iv);
% Algebraic expression distinguish positive from zero. [ae] is an
% AEX(c,c), [iv] an intervall containing the positive algebraic number
% represented by $ae$. Returns a RAT, which lies in $[0,ae]$
begin scalar rb;
rb := rat_mult(iv_rb iv,rat_mk(1,2));
%%% repeat ... until saves a sign test
while (aex_sgn aex_addrat(ae,rat_neg rb) neq 1) do
rb := rat_mult(rb,rat_mk(1,2));
return rb
end;
procedure aex_distinguishfromzero(ae,iv);
% Algebraic expression distinguish from zero. [ae] is an AEX(c,c),
% [iv] an intervall containing the algebraic number represented by
% $ae$. Returns a RAT, which lies in $[0,ae]$ or $[ae,0]$,
% respectively.
begin scalar sgnae,r;
sgnae := aex_sgn ae;
if eqn(sgnae,0) then
return rat_0();
r := if eqn(sgnae,1) then
aex_distinguishpositivefromzero(ae,iv)
else
aex_distinguishpositivefromzero(aex_neg ae,iv_neg iv);
return if eqn(sgnae,1) then r else rat_neg r
end;
procedure aex_cauchybound(ae,x);
% Algebraic expression cauchy bounds. [ae] is an univariate alg.
% polynomial in [x] with non-trivial leading coefficient. Returns
% an non-negative integer, the minimum of the cauchy bounds of ae.
begin scalar cfl,am,ctam,nb,ctl,ml,m,n,minabsam,cb,aesc;
if !*rlanuexdebug and aex_simplenullp aex_lc(ae,x) then
prin2t "***** aex_cauchybound: argument has trivial lc";
if aex_deg(ae,x) <= 0 then return rat_1(); % avoids trivial sturmchains
cfl := aex_coefl(ae,x); % has at least length 1
am := car cfl;
ctam := aex_containment am;
if iv_containszero ctam then <<
if ofsf_anuexverbosep() then
prin2 "+++ aex_cauchybound: iteration case +++";
nb := aex_distinguishfromzero(am,ctam);
ctam := if rat_less(nb,rat_0()) then
iv_mk(iv_lb ctam,nb)
else
iv_mk(nb,iv_rb ctam)
>>;
% now ctam is a containing interval for a_m without 0.
ctl := for each cf in cdr cfl collect aex_containment cf;
ml := for each iv in ctl collect iv_maxabs iv;
minabsam := iv_minabs ctam;
m := rat_max(rat_1(),rat_quot(rat_mapadd ml,minabsam));
n := if null ml then
rat_1()
else
rat_add(rat_1(),
rat_mapmax for each a in ml collect rat_quot(a,minabsam));
cb := rat_min(m,n);
aesc := aex_stdsturmchain(ae,x);
if !*rlanuexdebug and not (aex_sturmchainsgnch(aesc,x,cb) eq
aex_stchsgnch(aesc,x,'infty)) then
prin2 "aex_cauchybound: problem at right border";
if !*rlanuexdebug and not (aex_stchsgnch(aesc,x,'minfty) eq
aex_sturmchainsgnch(aesc,x,rat_minus(rat_neg cb,rat_1()))) then
prin2 "aex_cauchybound: problem at left border";
return cb;
end;
procedure aex_findrootsoflist(ael,x);
aex_findroots(aex_foldmult ael,x);
procedure aex_findroots(ae,x);
% Algebraic expression find roots. [ae] is an AEX(c,c+1). Returns a
% list of ANU(c+1). If [ae]'s ldeg is not positive, the empty list
% will be returned. The interval to start with has to be slightly
% enlarged.
begin scalar cb,rootlist;
if aex_deg(ae,x) < 1 then return nil;
%cb := rat_add(aex_cauchybound ae,rat_1()); % necessary to add 1.
cb := 256 . 1; %%% later cauchybound!!!
rootlist := aex_findrootsiniv1(ae,x,iv_mk(rat_neg cb,cb),
aex_stdsturmchain(ae,x));
return rootlist
end;
procedure aex_findrootsiniv(ae,x,iv);
% Algebraic expression find roots in interval. [ae] is an
% AEX(c,c+1), [iv] is an INT. The ldeg of ae has to be positive,
% i.e. [ae] must not represent a constant polynomial. Returns a
% ordered list of ANU(c+1).
aex_findrootsiniv1(ae,x,iv,aex_stdsturmchain(ae,x));
procedure aex_findrootsiniv1(ae,x,iv,sc);
% Algebraic expression find roots in interval. [ae] is an
% AEX(c,c+1), [iv] is an INT, [sc] is ae's sturmchain. Returns a
% ordered list of ANU(c+1). In particular, the ldeg of ae is
% positive, i.e. [ae].
begin scalar lb,rb,sclb,scrb,m,ml,mr,retl,r;
lb := iv_lb iv;
rb := iv_rb iv;
sclb := aex_sturmchainsgnch(sc,x,lb);
if sclb = 0 then return nil;
scrb := aex_sturmchainsgnch(sc,x,rb);
if sclb - scrb = 0 then
return nil;
if sclb - scrb = 1 then
return {anu_mk(ae,iv)};
m := rat_quot(rat_add(lb,rb),rat_fromnum 2);
ml := mr := m;
if aex_sgn aex_subrat1(ae,x,m) eq 0 then <<
r := aex_isoroot(ae,x,m,rat_mult(rat_minus(rb,lb),rat_mk(1,4)),sc);
ml := rat_minus(m,r);
mr := rat_add(m,r)
>>;
retl := aex_findrootsiniv1(ae,x,iv_mk(mr,rb),sc);
if not rat_eq(ml,mr) then
retl := anu_mk(ae,iv_mk(ml,mr)) . retl;
retl := append(aex_findrootsiniv1(ae,x,iv_mk(lb,ml),sc),retl);% nconc!!!
return retl
end;
procedure aex_isoroot(ae,x,m,r,sc);
% Algebraic expression isolate root. [ae] is an AEX(c,c+1), [m],
% [r] are RAT, [sc] is [ae]'s sturmchain. Returns a RAT $s$ such
% that [ae] has no root within the Interval $[m-s,m+s]$. In
% particular, the ldeg of ae is positive, i.e. [ae].
<< while not (aex_atrat0p(ae,x,rat_minus(m,r)) and
aex_atrat0p(ae,x,rat_add(m,r)) and
aex_sturmchainsgnch(sc,x,rat_minus(m,r))-
aex_sturmchainsgnch(sc,x,rat_add(m,r)) eq 1) do
r := rat_mult(r,rat_mk(1,2));
r
>>;
procedure aex_atrat0p(ae,x,r);
% Zero at rational point predicate. Check if [ae] is zero at [x].
aex_sgn aex_subrat1(ae,x,r) neq 0; %%% eq 0?
%%% --- global root isolation --- %%%
%procedure aex_globalrootisolationinit(ff,x);
% % [ff] is a (maybe empty) list of AEX (univariate alg.
% % polynomials). [ff] is factorized, i.e. each root occures only in
% % one polynomial.
% {nil,nil,for each f in ff collect f . aex_stdsturmchain(f,x)};
procedure aex_nextroot(rip,x);
% [rip] is of form {rootl,ivl,pscl}, where rootl is a list of ANU,
% ivl a list of IV, pscl a list of dotted pais p . sc where p ia an
% AEX (a univariate alg. polynomial) and sc is the corresponding
% sturmchain. Returns t, if a root is found, nil otherwise. The
% argument is changed in-place.
begin scalar rootfound,cb;
while not rootfound and rip_pscl rip do <<
% pscl is not empty. if ivl is empty, mk new iv with cb.
while null rip_ivl rip and rip_pscl rip do <<
cb := rat_add(aex_cauchybound(caar rip_pscl rip,x),rat_1());
rip_addivl(rip,iv_minuslist(iv_mk(rat_neg cb,cb),
for each a in rip_rootl rip collect anu_iv a));
% maybe after minuslist there are no intervals left
if null rip_ivl rip then rip_poppscl rip;
>>;
% care for the case that there was no poly left with a root
if rip_pscl rip then << % there is at least one interval and one poly
if !*rlanuexdebug and (null rip_ivl rip or null rip_pscl rip) then
prin2t "***** aex_nextroot: no interval or no polynomial";
if aex_nextroot1(rip,x) then rootfound := t;
if null rip_ivl rip then rip_poppscl rip
>>
>>;
if !*rlanuexdebug and rootfound then anu_check car rip_rootl rip;
if rootfound then return car rip_rootl rip else return nil
end;
procedure aex_nextroot1(rip,x);
% There is an interval and there is a poly. this function attempts
% to isolate at most one root. A root might be added to rootl and
% intervals might be added to ivl, but pscl is not popped. Returns
% t, if a root was found. Possible Optimazions: tagged intervals,
% use anu_fromrat if a rat. root is found
begin scalar iv,lb,rb,sc,sclb,scrb,m,ml,mr,r;
% remove iv from ivl
iv := rip_popivl rip; lb := iv_lb iv; rb := iv_rb iv;
sc := cdar rip_pscl rip;
sclb := aex_sturmchainsgnch(sc,x,lb);
if sclb = 0 then return nil;
scrb := aex_sturmchainsgnch(sc,x,rb);
if sclb - scrb = 0 then return nil;
if sclb - scrb = 1 then <<
rip_addroot(rip,aex_refinewrt1(anu_mk(caar rip_pscl rip,iv),
cdr rip_pscl rip,cdar rip_pscl rip));
return t
>>;
% there are at least two roots.
m := rat_quot(rat_add(lb,rb),rat_fromnum 2);
ml := mr := m;
if aex_sgn aex_subrat1(caar rip_pscl rip,x,m) eq 0 then <<
r := aex_isoratroot(m,rat_mult(rat_minus(rb,lb),rat_mk(1,4)),
rip_pscl rip,x); % cdr rip_pscl rip would be wrong
ml := rat_minus(m,r); mr := rat_add(m,r);
%rip_addroot(rip,anu_mk(caar rip_pscl rip,iv_mk(ml,mr)))
% instead, the following is only a very small optimization(1-2%):
rip_addroot(rip,anu_fromrat(x,m,iv_mk(ml,mr)))
>>;
rip_addivl(rip,{iv_mk(mr,rb)});
rip_addivl(rip,{iv_mk(lb,ml)});
return not rat_eq(ml,mr);
end;
procedure aex_isoratroot(m,r,pscl,x);
% [m] is a rational root of [caar pscl]. Returns a r such that no
% sturmchain from cdr scl has a sign change in (m-r,m+r), and such
% that car scl has 1 sign change in (m-r,m+r).
<< % Refine r until p's sc has only 1 sgn change.
while aex_atratnullp(caar pscl,x,rat_minus(m,r)) or
aex_atratnullp(caar pscl,x,rat_minus(m,r)) or
aex_deltastchsgnch(cdar pscl,x,
iv_mk(rat_minus(m,r),rat_add(m,r))) > 1 do
r := rat_mult(r,rat_mk(1,2));
pscl := cdr pscl;
% Refine r until the other sc's have no sgn changes.
for each psc in pscl do
while aex_atratnullp(car psc,x,rat_minus(m,r)) or
aex_atratnullp(car psc,x,rat_minus(m,r)) or
aex_deltastchsgnch(cdr psc,x,
iv_mk(rat_minus(m,r),rat_add(m,r))) > 0 do
r := rat_mult(r,rat_mk(1,2));
r
>>;
procedure aex_refinewrt1(alpha,pscl,scalpha);
% Refine [alpha] st. the sturm chain [scl] has no sign change in
% alpha's isolating interval. scalpha is the sturm chain for
% alpha's defining polynomial.
<< for each psc in pscl do %%% null test necessary
while aex_atratnullp(car psc,x,iv_lb anu_iv alpha) or
aex_atratnullp(car psc,x,iv_rb anu_iv alpha) or
aex_deltastchsgnch(cdr psc,x,anu_iv alpha) > 0 do
anu_refine1ip(alpha,scalpha);
alpha
>> where x=aex_mvar anu_dp alpha;
procedure aex_deltastchsgnch(sc,x,iv);
% delta sturm chain sign changes in an interval
begin scalar sclb,scrb;
sclb := aex_sturmchainsgnch(sc,x,iv_lb iv);
if sclb eq 0 then return 0;
scrb := aex_sturmchainsgnch(sc,x,iv_rb iv);
return sclb - scrb;
end;
procedure aex_atratnullp(ae,x,r);
aex_sgn aex_subrat1(ae,x,r) eq 0;
% - little helpers for global root isolation - %
% rip stands for root list, interval list, p.sc list
procedure psc_sortdesc(psc1,psc2);
% Sort function descending. psc1 and psc2 are list of dotted pairs
% p.sc phere p is a ratpoly and sc is a sturmchain. To be used by
% generic_sort, yiels a sorted list with descending length of sturm
% chains.
sgn(length cdr psc2 - length cdr psc1);
procedure psc_sortasc(psc1,psc2);
% Sort function asscending. psc1 and psc2 are list of dotted pairs
% p.sc phere p is a ratpoly and sc is a sturmchain. To be used by
% generic_sort, yiels a sorted list with ascending length of sturm
% chains.
sgn(length cdr psc1 - length cdr psc2);
procedure rip_init(ael,x);
% [ael] is a (maybe empty) list of AEX (univariate alg.
% polynomials). [ael] is factorized, i.e. each root occures only in
% one polynomial.
{nil,nil,for each ae in ael collect ae . aex_stdsturmchain(ae,x)};
% generic_sort(for each ae in ael collect ae . aex_stdsturmchain(ae,x),
% 'psc_sortasc)};
procedure rip_rootl rip;
car rip;
procedure rip_ivl rip;
cadr rip;
procedure rip_pscl rip;
caddr rip;
procedure rip_putivl(rip,ivl);
cadr rip := ivl;
procedure rip_putpscl(rip,pscl);
caddr rip := pscl;
procedure rip_addroot(rip,root);
car rip := root . car rip;
procedure rip_addivl(rip,ivl);
cadr rip := nconc(ivl,rip_ivl rip);
procedure rip_poppscl rip;
caddr rip := cdr caddr rip;
procedure rip_popivl rip;
<<cadr rip := cdr cadr rip; carivl>> where carivl=car rip_ivl rip;
%procedure aex_ivlfromrootl rootl;
% % interval list from root list. [rootl] is a list of ANU.
% for each a in rootl collect anu_iv a;
% module anu;
procedure anu_mk(ae,iv);
{'anu,ae,iv};
procedure anu_dp(a);
cadr a;
procedure anu_iv(a);
caddr a;
procedure anu_putiv(a,iv);
% Algebraic number get interval. [a] is an ANU. [iv] is an
% INT. Changes the isolating interval in a nonconstructively.
% The returned value is not of interest.
caddr a := iv;
procedure anu_print a;
<<
prin2 "("; aex_print anu_dp a; prin2 ", ";
iv_print anu_iv a; prin2 ")";
>>;
procedure anu_check(a);
% Algebraic number check.
begin scalar dp,x,l,r,s,valid;
valid := t;
dp := anu_dp a;
x := aex_freeids dp; % tmp
if length x neq 1 then
prin2t "***** anu_check: def. poly corrupt";
x := car x; l := iv_lb anu_iv a; r := iv_rb anu_iv a;
s := aex_stdsturmchain(dp,aex_mvar dp);
if aex_nullp aex_subrat1(dp,x,l) then <<
valid := nil;
prin2t "***** anu_check: def. poly null at left bound";
>>;
if aex_nullp aex_subrat1(dp,x,r) then <<
valid := nil;
prin2t "***** anu_check: def. poly null at right bound";
>>;
if aex_deltastchsgnch(s,x,anu_iv a) neq 1 then <<
valid := nil;
prin2t "***** anu_check: no root";
>>;
return valid;
end;
procedure anu_fromrat(x,r,iv);
anu_mk(aex_fromrp ratpoly_mklin(x,r),iv);
procedure anu_refine1ip(a,s);
% Algebraic number refine in place . [a] is an ANU, [s] is a
% sturmchain. Returns an [a] with smaller isolating intervall.
% Remark: 1000 replaced by 4.
begin scalar x,iv,lb,rb,m,scm;
x := aex_mvar anu_dp a;
iv := anu_iv a;
lb := iv_lb iv; rb := iv_rb iv;
m := rat_quot(rat_add(lb,rb),rat_mk(2,1));
if aex_sgn aex_subrat1(anu_dp a,x,m) = 0 then <<
anu_putiv(a,iv_mk(
rat_minus(m,rat_mult(rat_mk(1,4),rat_minus(m,lb))),
rat_add(m,rat_mult(rat_mk(1,4),rat_minus(m,lb)))));
return a;
>>;
scm := aex_sturmchainsgnch(s,x,m);
if (aex_sturmchainsgnch(s,x,lb)-scm) = 1 then <<
anu_putiv(a,iv_mk(lb,m));
return a;
>>;
anu_putiv(a,iv_mk(m,rb));
return a
end;
procedure sf_idl f;
% if there is a main variable, then it will be the car.
if not domainp f then
mvar f . setunion(sf_idl lc f,setminus(sf_idl red f,{mvar f}));
procedure setminus(ss1,ss2);
for each s1 in ss1 join if not member(s1,ss2) then {s1};
procedure intersection(ss1,ss2);
setminus(ss1,setminus(ss1,ss2));
%procedure setunion(ss1,ss2);
% append(setminus(ss1,ss2),ss2);
procedure idorder();
begin scalar idorder;
idorder := setkorder nil;
setkorder idorder;
return idorder;
end;
endmodule; % ofsfanuex
end; % of file