login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit
Author
Mike Dewar
Date Created
November 1996
Description:
This category is an implementation of interval arithmetic and transcendental functions over intervals.
Note:
It appears to be nearly identical to the SPAD category interval.spad.

aldor
#include "axiom.as"
FUNCAT ==> Join(FloatingPointSystem,TranscendentalFunctionCategory);
define IntervalCategory2(R:FUNCAT): Category ==
RetractableTo(Integer))
with {
approximate;
interval : (R,R) -> %;
++ interval(inf,sup) creates a new interval, either \axiom{[inf,sup]} if
++ \axiom{inf <= sup} or \axiom{[sup,in]} otherwise.
qinterval : (R,R) -> %;
++ qinterval(inf,sup) creates a new interval \axiom{[inf,sup]}, without
++ checking the ordering on the elements.
interval : R -> %;
++ interval(f) creates a new interval around f.
interval : Fraction Integer -> %;
++ interval(f) creates a new interval around f.
inf : % -> R;
++ inf(u) returns the infinum of \axiom{u}.
sup : % -> R;
++ sup(u) returns the supremum of \axiom{u}.
width : % -> R;
++ width(u) returns \axiom{sup(u) - inf(u)}.
positive? : % -> Boolean;
++ positive?(u) returns \axiom{true} if every element of u is positive,
++ \axiom{false} otherwise.
negative? : % -> Boolean;
++ negative?(u) returns \axiom{true} if every element of u is negative,
++ \axiom{false} otherwise.
contains? : (%,R) -> Boolean;
++ contains?(i,f) returns true if \axiom{f} is contained within the interval
++ \axiom{i}, false otherwise.
}
import from Integer;
import from R;
Rep ==> Record(Inf:R, Sup:R);
import from Rep;
local roundDown(u:R):R ==
if zero?(u) then float(-1,-(bits() pretend Integer));
else float(mantissa(u) - 1,exponent(u));
local roundUp(u:R):R   ==
if zero?(u) then float(1, -(bits()) pretend Integer);
else float(mantissa(u) + 1,exponent(u));
-- Sometimes the float representation does not use all the bits (e.g. when
-- representing an integer in software using arbitrary-length Integers as
-- your mantissa it is convenient to keep them exact).  This function
-- normalises things so that rounding etc. works as expected.  It is only
-- called when creating new intervals.
local normaliseFloat(u:R):R ==
if zero? u then u else {
m : Integer := mantissa u;
b : Integer := bits() pretend Integer;
l : Integer := length(m);
if (l < b) then {
BASE : Integer := base()$R pretend Integer; float(m*BASE**((b-l) pretend PositiveInteger),exponent(u)-b+l); } else u; } interval(i:R,s:R):% == { i > s => per [roundDown normaliseFloat s,roundUp normaliseFloat i]; per [roundDown normaliseFloat i,roundUp normaliseFloat s]; } interval(f:R):% == { zero?(f) => 0; one?(f) => 1; -- This next part is necessary to allow e.g. mapping between Expressions: -- AXIOM assumes that Integers stay as Integers! import from Union(value1:Integer,failed:'failed'); fnew : R := normaliseFloat f; retractIfCan(f)@Union(value1:Integer,failed:'failed') case value1 => per [fnew,fnew]; per [roundDown fnew, roundUp fnew]; } qinterval(i:R,s:R):% == per [roundDown normaliseFloat i,roundUp normaliseFloat s]; local exactInterval(i:R,s:R):% == per [i,s]; local exactSupInterval(i:R,s:R):% == per [roundDown i,s]; local exactInfInterval(i:R,s:R):% == per [i,roundUp s]; inf(u:%):R == (rep u).Inf; sup(u:%):R == (rep u).Sup; width(u:%):R == (rep u).Sup - (rep u).Inf; contains?(u:%,f:R):Boolean == (f > inf(u)) and (f < sup(u)); positive?(u:%):Boolean == inf(u) > 0; negative?(u:%):Boolean == sup(u) < 0; (<)(a:%,b:%):Boolean == if inf(a) < inf(b) then true else if inf(a) > inf(b) then false else sup(a) < sup(b); (+)(a:%,b:%):% == { -- A couple of blatent hacks to preserve the Ring Axioms! if zero?(a) then return(b) else if zero?(b) then return(a); if a=b then return qinterval(2*inf(a),2*sup(a)); qinterval(inf(a) + inf(b), sup(a) + sup(b)); } (-)(a:%,b:%):% == { if zero?(a) then return(-b) else if zero?(b) then return(a); if a=b then 0 else qinterval(inf(a) - sup(b), sup(a) - inf(b)); } (*)(a:%,b:%):% == { -- A couple of blatent hacks to preserve the Ring Axioms! if one?(a) then return(b) else if one?(b) then return(a); if zero?(a) then return(0) else if zero?(b) then return(0); prods : List R := sort [inf(a)*inf(b),sup(a)*sup(b), inf(a)*sup(b),sup(a)*inf(b)]; qinterval(first prods, last prods); } (*)(a:Integer,b:%):% == { if (a > 0) then qinterval(a*inf(b),a*sup(b)); else if (a < 0) then qinterval(a*sup(b),a*inf(b)); else 0; } (*)(a:PositiveInteger,b:%):% == qinterval(a*inf(b),a*sup(b)); (**)(a:%,n:PositiveInteger):% == { contains?(a,0) and zero?((n pretend Integer) rem 2) => interval(0,max(inf(a)**n,sup(a)**n)); interval(inf(a)**n,sup(a)**n); } (^) (a:%,n:PositiveInteger):% == { contains?(a,0) and zero?((n pretend Integer) rem 2) => interval(0,max(inf(a)**n,sup(a)**n)); interval(inf(a)**n,sup(a)**n); } (-)(a:%):% == exactInterval(-sup(a),-inf(a)); (=)(a:%,b:%):Boolean == (inf(a)=inf(b)) and (sup(a)=sup(b)); (~=)(a:%,b:%):Boolean == (inf(a)~=inf(b)) or (sup(a)~=sup(b)); 1:% == {one : R := normaliseFloat 1; per([one,one])}; 0:% == per([0,0]); recip(u:%):Union(value1:%,failed:'failed') == { contains?(u,0) => [failed]; vals:List R := sort[1/inf(u),1/sup(u)]; [qinterval(first vals, last vals)]; } unit?(u:%):Boolean == contains?(u,0); exquo(u:%,v:%):Union(value1:%,failed:'failed') == { contains?(v,0) => [failed]; one?(v) => [u]; u=v => [1]; u=-v => [-1]; vals:List R := sort[inf(u)/inf(v),inf(u)/sup(v),sup(u)/inf(v),sup(u)/sup(v)]; [qinterval(first vals, last vals)]; } gcd(u:%,v:%):% == 1; coerce(u:Integer):% == { ur := normaliseFloat(u::R); exactInterval(ur,ur); } interval(u:Fraction Integer):% == { import { log2 : % -> %; coerce : Integer -> %; retractIfCan : % -> Union(value1:Integer,failed:'failed');} from Float; flt := u::R; -- Test if the representation in R is exact --den := denom(u)::Float; local bin : Union(value1:Integer,failed:'failed'); bin := retractIfCan(log2(denom(u)::Float)); bin case value1 and length(numer u)$Integer < (bits() pretend Integer) => {
flt := normaliseFloat flt;
exactInterval(flt,flt);
}
qinterval(flt,flt);
}
retractIfCan(u:%):Union(value1:Integer,failed:'failed') == {
not zero? width(u) => [failed];
retractIfCan inf u;
}
retract(u:%):Integer == {
not zero? width(u) =>
error "attempt to retract a non-Integer interval to an Integer";
retract inf u;
}
coerce(u:%):OutputForm ==
bracket([coerce inf(u), coerce sup(u)]$List(OutputForm)); characteristic():NonNegativeInteger == 0; -- Explicit export from TranscendentalFunctionCategory pi():% == qinterval(pi(),pi()); -- From ElementaryFunctionCategory log(u:%):% == { positive?(u) => qinterval(log inf u, log sup u); error "negative logs in interval"; } exp(u:%):% == qinterval(exp inf u, exp sup u); (**)(u:%,v:%):% == { zero?(v) => if zero?(u) then error "0**0 is undefined" else 1; one?(u) => 1; expts : List R := sort [inf(u)**inf(v),sup(u)**sup(v), inf(u)**sup(v),sup(u)**inf(v)]; qinterval(first expts, last expts); } -- From TrigonometricFunctionCategory -- This function checks whether an interval contains a value of the form -- offset + 2 n pi'. local hasTwoPiMultiple(offset:R,Pi:R,i:%):Boolean == { import from Integer; next : Integer := retract ceiling( (inf(i) - offset)/(2*Pi) ); contains?(i,offset+2*next*Pi); } -- This function checks whether an interval contains a value of the form -- offset + n pi'. local hasPiMultiple(offset:R,Pi:R,i:%):Boolean == { import from Integer; next : Integer := retract ceiling( (inf(i) - offset)/Pi ); contains?(i,offset+next*Pi); } sin(u:%):% == { import from Integer; Pi : R := pi(); hasOne? : Boolean := hasTwoPiMultiple(Pi/(2::R),Pi,u); hasMinusOne? : Boolean := hasTwoPiMultiple(3*Pi/(2::R),Pi,u); if hasOne? and hasMinusOne? then exactInterval(-1,1); else { vals : List R := sort [sin inf u, sin sup u]; if hasOne? then exactSupInterval(first vals, 1); else if hasMinusOne? then exactInfInterval(-1,last vals); else qinterval(first vals, last vals); } } cos(u:%):% == { Pi : R := pi(); hasOne? : Boolean := hasTwoPiMultiple(0,Pi,u); hasMinusOne? : Boolean := hasTwoPiMultiple(Pi,Pi,u); if hasOne? and hasMinusOne? then exactInterval(-1,1); else { vals : List R := sort [cos inf u, cos sup u]; if hasOne? then exactSupInterval(first vals, 1); else if hasMinusOne? then exactInfInterval(-1,last vals); else qinterval(first vals, last vals); } } tan(u:%):% == { Pi : R := pi(); if width(u) > Pi then error "Interval contains a singularity" else { -- Since we know the interval is less than pi wide, monotonicity implies -- that there is no singularity. If there is a singularity on a endpoint -- of the interval the user will see the error generated by R. lo : R := tan inf u; hi : R := tan sup u; lo > hi => error "Interval contains a singularity"; qinterval(lo,hi); } } csc(u:%):% == { Pi : R := pi(); if width(u) > Pi then error "Interval contains a singularity" else { import from Integer; -- singularities are at multiples of Pi if hasPiMultiple(0,Pi,u) then error "Interval contains a singularity"; vals : List R := sort [csc inf u, csc sup u]; if hasTwoPiMultiple(Pi/(2::R),Pi,u) then exactInfInterval(1,last vals); else if hasTwoPiMultiple(3*Pi/(2::R),Pi,u) then exactSupInterval(first vals,-1); else qinterval(first vals, last vals); } } sec(u:%):% == { Pi : R := pi(); if width(u) > Pi then error "Interval contains a singularity" else { import from Integer; -- singularities are at Pi/2 + n Pi if hasPiMultiple(Pi/(2::R),Pi,u) then error "Interval contains a singularity"; vals : List R := sort [sec inf u, sec sup u]; if hasTwoPiMultiple(0,Pi,u) then exactInfInterval(1,last vals); else if hasTwoPiMultiple(Pi,Pi,u) then exactSupInterval(first vals,-1); else qinterval(first vals, last vals); } } cot(u:%):% == { Pi : R := pi(); if width(u) > Pi then error "Interval contains a singularity" else { -- Since we know the interval is less than pi wide, monotonicity implies -- that there is no singularity. If there is a singularity on a endpoint -- of the interval the user will see the error generated by R. hi : R := cot inf u; lo : R := cot sup u; lo > hi => error "Interval contains a singularity"; qinterval(lo,hi); } } -- From ArcTrigonometricFunctionCategory asin(u:%):% == { lo : R := inf(u); hi : R := sup(u); if (lo < -1) or (hi > 1) then error "asin only defined on the region -1..1"; qinterval(asin lo,asin hi); } acos(u:%):% == { lo : R := inf(u); hi : R := sup(u); if (lo < -1) or (hi > 1) then error "acos only defined on the region -1..1"; qinterval(acos hi,acos lo); } atan(u:%):% == qinterval(atan inf u, atan sup u); acot(u:%):% == qinterval(acot sup u, acot inf u); acsc(u:%):% == { lo : R := inf(u); hi : R := sup(u); if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then error "acsc not defined on the region -1..1"; qinterval(acsc hi, acsc lo); } asec(u:%):% == { lo : R := inf(u); hi : R := sup(u); if ((lo < -1) and (hi > -1)) or ((lo < 1) and (hi > 1)) then error "asec not defined on the region -1..1"; qinterval(asec lo, asec hi); } -- From HyperbolicFunctionCategory tanh(u:%):% == qinterval(tanh inf u, tanh sup u); sinh(u:%):% == qinterval(sinh inf u, sinh sup u); sech(u:%):% == { negative? u => qinterval(sech inf u, sech sup u); positive? u => qinterval(sech sup u, sech inf u); vals : List R := sort [sech inf u, sech sup u]; exactSupInterval(first vals,1); } cosh(u:%):% == { negative? u => qinterval(cosh sup u, cosh inf u); positive? u => qinterval(cosh inf u, cosh sup u); vals : List R := sort [cosh inf u, cosh sup u]; exactInfInterval(1,last vals); } csch(u:%):% == { contains?(u,0) => error "csch: singularity at zero"; qinterval(csch sup u, csch inf u); } coth(u:%):% == { contains?(u,0) => error "coth: singularity at zero"; qinterval(coth sup u, coth inf u); } -- From ArcHyperbolicFunctionCategory acosh(u:%):% == { inf(u)<1 => error "invalid argument: acosh only defined on the region 1.."; qinterval(acosh inf u, acosh sup u); } acoth(u:%):% == { lo : R := inf(u); hi : R := sup(u); if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then error "acoth not defined on the region -1..1"; qinterval(acoth hi, acoth lo); } acsch(u:%):% == { contains?(u,0) => error "acsch: singularity at zero"; qinterval(acsch sup u, acsch inf u); } asech(u:%):% == { lo : R := inf(u); hi : R := sup(u); if (lo <= 0) or (hi > 1) then error "asech only defined on the region 0 < x <= 1"; qinterval(asech hi, asech lo); } asinh(u:%):% == qinterval(asinh inf u, asinh sup u); atanh(u:%):% == { lo : R := inf(u); hi : R := sup(u); if (lo <= -1) or (hi >= 1) then error "atanh only defined on the region -1 < x < 1"; qinterval(atanh lo, atanh hi); } -- From RadicalCategory (**)(u:%,n:Fraction Integer):% == interval(inf(u)**n,sup(u)**n); } aldor  Compiling FriCAS source code from file /var/zope2/var/LatexWiki/interval2.as using AXIOM-XL compiler and options -O -Fasy -Fao -Flsp -laxiom -Mno-ALDOR_W_WillObsolete -DAxiom -Y$AXIOM/algebra -I $AXIOM/algebra Use the system command )set compiler args to change these options. "/var/zope2/var/LatexWiki/interval2.as", line 9: approximate; ..^ [L9 C3] #1 (Error) No meaning for identifier approximate'. "/var/zope2/var/LatexWiki/interval2.as", line 66: float(m*BASE**((b-l) pretend PositiveInteger),exponent(u)-b+l); ..............^ [L66 C15] #2 (Error) Argument 1 of **' did not match any possible parameter type. The rejected type is Integer. Expected type %. "/var/zope2/var/LatexWiki/interval2.as", line 147: interval(0,max(inf(a)**n,sup(a)**n)); .....................^.........^ [L147 C22] #3 (Error) Argument 1 of **' did not match any possible parameter type. The rejected type is R. Expected type %. "/var/zope2/var/LatexWiki/interval2.as", line 148: interval(inf(a)**n,sup(a)**n); .............^.........^ [L148 C14] #5 (Error) Argument 1 of **' did not match any possible parameter type. The rejected type is R. Expected type %. "/var/zope2/var/LatexWiki/interval2.as", line 153: interval(0,max(inf(a)**n,sup(a)**n)); .....................^.........^ [L153 C22] #7 (Error) Argument 1 of **' did not match any possible parameter type. The rejected type is R. Expected type %. "/var/zope2/var/LatexWiki/interval2.as", line 154: interval(inf(a)**n,sup(a)**n); .............^.........^ [L154 C14] #9 (Error) Argument 1 of **' did not match any possible parameter type. The rejected type is R. Expected type %. [L154 C24] #11 (Fatal Error) Too many errors (use -M emax=n' or -M no-emax' to change the limit). The )library system command was not called after compilation. axiom x:=interval(1.1,2.2)$Interval2(Float)
There are no library operations named Interval2
Use HyperDoc Browse or issue
)what op Interval2
to learn if there is any operation containing " Interval2 " in
its name.
Cannot find a definition or applicable library operation named
Interval2 with argument type(s)
Domain
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need. sin(x)  (1) Type: Expression(Integer) Reason that it looks like interval.spad is that I rewrote interval.as for the SPAD compiler, so that we are able to run Axiom without Aldor. Juergen Weiss Juergen Weiss wrote: I rewrote interval.as for the SPAD compiler Oh yes, I am sorry that I forgot. I think your important contribution should be so noted in the interval.spad.pamphlet file! Here's a patch: $ diff -Nar interval.spad.pamphlet.orig interval.spad.pamphlet
5c5
< \author{Mike Dewar}
---
> \author{Juergen Weiss}
7a8,9
> This is a re-write of interval.as' in SPAD to enable all of
> Axiom's library to be compiled without the Aldor compiler.
15c17
< +++ Author: Mike Dewar
---
> +++ Original Author: Mike Dewar
`

--------

 Subject:   Be Bold !! ( 14 subscribers )