\documentclass{article}
\usepackage{fullpage}
\title{Interval Arithmetic for Maxima: A Brief Summary}
\author{Richard J. Fateman\\
Computer Science Division\\
Electrical Engineering and Computer Sciences\\
University of California at Berkeley }
\begin{document}
\maketitle
\begin{abstract}
%May, June, July, 2006
%March, 2013
%Sept 2013 touchups mostly
%June 2014 added notes about Maxima and property lists.
%Nov 2015, revived CAS as the way to do stuff toward fixing the dependency
%problems, SUE production etc. And unknown other patches to make
%interval arithmetic work better.
%March 2016, cut it all out except the most basic stuff
This is a short introduction to an implementation
of interval arithmetic data structures and operations
for the Maxima computer algebra
system.
We expect that further algorithms necessary to implement the
forthcoming IEEE 1788 standard for interval arithmetic can be readily
fit into this structure.
\end{abstract}
\section {Introduction}
Computer algebra systems (CAS) have a quite flexible way of introducing
new concepts into their mathematical systems: functional notation.
In particular, it is tempting to just introduce an interval which
mathematicians might write as [a,b] in Maxima
by typing {\tt interval(a,b)}.
Using this most obvious representation of an interval
and using the usual techniques for programming operations with it results in some
uncomfortable consequences.
In particular, one encounters the so-called dependency problem.
We can fix this and then proceed to build one (or several different)
models for intervals that are more convenient for users.
\section{Why the obvious method is unsatisfactory}
If {\tt interval(1,2)} is an implementation
for the interval consisting of some $x$ such that $1 \leq x \leq 2$.
then {\tt interval(1,2) - interval(1,2)} is,
by the standard simplification for {\tt x - x}, reduced to 0. This is correct
if the computation is of $s-s$ where each $s$ is {\em the same} {\tt interval(1,2)}.
On the other hand, it is not zero but {\tt interval(-1,1)}
if the computation is of $s-t$ where $t$ is some
{\em other}
$y$ such that $1 \leq y \leq 2$. That is, $t$ is independent of $s$.
These kinds of results affect other operations as well.
\section{A solution}
There are many solutions possible, see Fateman\cite{fate} for
a more expansive discussion of alternative representations in
an algebra system. Also, we may wish to support some of the
numerical tools which
can be found on internet for interval computations.
However, our context is different. We are not restricted to
endpoints consisting of machine-precision floating-point
approximations. We are capable of manipulating symbolic expressions
which are later committed to interval evaluation, e.g.
$9\,x^2+24\,x+9$ is inferior, as an interval expression,
to $\left(4+3\,x\right)^2-7$.
We would also like to be compatible with the forthcoming
IEEE standard 1788 for interval arithmetic\footnote{
http://grouper.ieee.org/groups/1788/}.
Our favorite solution to the dependence problem is to use Lisp
structures for intervals. The Maxima system sees the intervals, at
least ``at first glance'', as atomic symbols. Thus we can define a program
to be used as {\tt ri(1,2)} to constructs an ``atomic'' real interval
representing $x$ for $1\leq x \leq 2$. Uttering {\tt ri(1,2)} a
second time constructs a {\em different} atomic object that happens to
have the same values (1,2) stored internally. Maxima simplifies the difference of
two identical structures (that is, stored in the same memory location) to
zero, but not if two such intervals are independent. The independence
problem is solved.
Arithmetic is defined in the usual prescribed manner first by
using accessor functions to extract upper and lower limits from the
operands, perform arithmetic (with rounding as appropriate) and build
new interval objects. A method for printing must also be
supplied, and we have decided that it would be better to have
an {\em unambiguously different} display for intervals.
We have programmed {\tt <1..2>} to be the display format for the
interval produced by (for example) typing {\tt ri(1,2)}.
\section{Endpoints}
We are considering {\em real intervals}, and so the form {\tt
ri($a,~b$)}
is acceptable only if {$ a,~ b$} are Maxima (real) numbers, with
$a**b$
is empty; to be specific, 1788 suggests $[\infty,-\infty]$. A NaI may
involve floating-point NaNs (Not a Number) but the standard does not
constrain the format.
There is an important issue that arises, since a package to do
operations on intervals must do more that arithmetic on intervals {\em
per se}. Among these: we can try to extract an endpoint from an
``improper'' interval {\tt <-oo..oo>}. What are we to see? One
possibility is to signal an exception. This seems harsh at first
glance, but if we can test for improper intervals perhaps we can
side-step this possibility. Even more extreme, note that extracting
left/right endpoints from NaI cannot yield any regular number, so a
careful program must either check for the NaI (again, extracting
an endpoint results in an exception),
or later check for the value (whatever is extracted)
before the next step, which might be arithmetic or comparison.
{If we return one of the IEEE 754 floating-point number
``infinities'', or for NaI, a NaN, we then must support further operations
on them. Since intervals may not even involve floats, as for example
z=[-1/2, 1/3], why should 1/z involve floats? The endpoints of z are
rational, so should the endpoints of 1/z be rational infinities? Are
we talking nonsense here?}
Perhaps we should return symbols, which have been
used for decades in Maxima, but whose
support is spotty: ``inf'' or ``infinity''.
(Note: Maxima, in its standard configuration
will simplify ``inf-inf'' to zero.) Consequently, our
current design is to return an indexed
object such as {\tt infinity[n]} with a running, automatically-incremented
index.
{There are other choices we have examined and implemented
including changing the system to allow 1/0 as acceptable data to
return as representing
infinity. In such a system with related other rational-like objects,
we can also construct something akin to ``not a number'' by using
0/0. It would be trivial to make this change to the interval
data design, but there are ramifications for existing programs.}
In particular, this rational number style
frees the representation of $\infty$ from the domain of IEEE floating-point but requires
that we insert new checks in Maxima's simplifier code to
treat such ``rational infinities'' appropriately. For example
1/0-1/0 returns 0/0 etc.
Computing with infinities {\em outside intervals} in any of several
designs for Maxima is discussed in another paper, forthcoming.
\section{Next} Beyond the current level of implementation
we can incrementally build a whole host of
library routines for operating on intervals, using the IEEE 1788
standard as a template when appropriate.
Additionally, a computer algebra system can support quite different
operations from those of a numerical library. As one example,
programs striving
to rearrange computations as
``single use expressions'' (SUE) \cite{walster}, or other beneficial
forms, or increasing precision and/or subdividing intervals for
improved
results.
A sample version of an
interval trigonometric routine (sine) is provided, to exhibit an
approach using Lisp. Writing in Lisp is not at all required though.
Models of such routines from libraries written in other (typically
C-like)
languages can either be loaded into Lisp, or can be effectively
transliterated into
Lisp language programs. The latter case may be preferable
if the library does not implement arbitrary-precision arithmetic.
In the example below, {\tt f()} computes the negation of a proper
interval, K[a,b] by returning -K= [-b,-a]
\begin{verbatim}
matchdeclare(anint,ri_p)$
/* In a matching rule or "tellsimp" the variable anint matches an object that passes
the test ri_p */
tellsimp(f(anint),ri(-ri_hi(anint),-ri_lo(anint)));
/* for example, f(ri(1,3)), returns ri(-3,-1) which prints as <-3..-1>
\end{verbatim}
A more bulletproof version of this program would have to check for
NaI, and (perhaps) worry about infinite endpoints.
\section{How to get this into your system}
The program features are defined in this file:
\verb|http://www.cs.berkeley.edu/~fateman/lisp/interval.lisp|
which can be downloaded and loaded into
any
Maxima system. This defines functions
{\tt ri, ri\_p, ri\_lo or left, ri\_hi or right, widen\_ri, interval\_within,
interval\_intersection ?bu, ?bd, inf, minf, ind,
nai\_p, ri\_lo\_minf\_p, ri\_hi\_inf\_p, doint.}
Expressions like {\tt H: ri(3,4)-r(3,4)} are left unchanged, but the
{\tt doint} command ``does interval simplification''. Therefore
{\tt doint(H)} results in an {\tt <-1..1>}. Sometimes the
expression can simplified without {\tt doint}, as for example,
sin(interval). The rationale for {\tt doint} is that we do not
wish to burden the general simplifier with the task of incessantly
checking for intervals\footnote{An additional rationale is that
we do not wish to burden the author/programmer with the task of
re-programming so much of the Maxima general simplifier.}.
{\tt ri(a,b)}, where {\tt a} and {\tt b} are numeric creates a real
interval.
Note,
{\tt ri(1,2)+x} does {\em not} produce {\tt<1+x,2+x>} but is
left alone. As mentioned earlier, the utility of symbolic
endpoints for intervals is doubtful.
{\tt widen\_ri(a,b)} produces an interval wider than {\tt ri(a,b)} if
it is given floating-point numbers, by rounding outward.
Exact numbers are not widen. A single argument {\tt widen\_ri(a)}
bumps the single float down and up.
{?bu} and {?bd} for bump up and bump provide a facility
that is general but slightly cruder than we would
like. It would be technically preferable for machine floats -- that is, it would produce
tighter bounds, if we could wrap operations with
floating-point rounding modes as {\tt block([?round\_mode:up],
a+b) }
but this is not supportable using only Common Lisp standard
operations.
Consequently we perform the operations and then bump the result.
This loss in tightness (at the cost of speed) can be remedied to some extent
by using software bigfloats of higher precision for the endpoints.
The Lisp tradition is to indicate programs that return true or false,
(predicates) with a name ending in {\tt -p}. This use of ``-''
is inconvenient, so we use {\tt \_p} as in the predicate testing
to see if an object is a real interval, {\tt ri\_p}.
The other operations are obvious from their names.
\section{Further information}
The bibliography provides links to many more detailed
commentaries on interval arithmetic, as does the IEEE 1788
standard. This brief paper indicates a bare minimum of
information on the program code available in Maxima.
%\verb|q:ri(-2,3), doint(1/q) is <-oo..oo>. left(%)| is -infinity[8],
%some indexed infinity.
\begin{thebibliography}{99}
\bibitem {berz-taylor} Berz, M., COSY INFINITY system, http://bt.pa.msu.edu/index\_cosy.htm
%{\em Reliable Computing 4:} 83-97.
\bibitem{emiris} Ioannis Emiris, Richard Fateman,
``Towards an Efficient Implementation of Interval Arithmetic,''
Technical Report UCB/CSD-92-93, July, 1992 (UC Berkeley)
{\tt http://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6247.html}
\bibitem{fate} Richared Fateman, ``To Infinity and Beyond'',
{\verb|http://www.cs.berkeley.edu/~fateman/papers/infinity.pdf|}
\bibitem{hansen}
Eldon R. Hansen, ``Sharpness in Interval Computations'',
{\em Reliable Computing} 3: 17-29 (1997).
\bibitem{p1788} IEEE,
{\tt http://standards.ieee.org/findstds/standard/1788-2015.html}
\bibitem {martin}
Ralph Martin, Huahao Shou, Irina Voiculescu, Adrian Bowyer and Guojin
Wang, ``Comparison of interval methods for plotting algebraic curves,''
Computer Aided Geometric Design, Volume 19, Issue 7, July 2002,
Pages 553-587.
\bibitem{generic} Richard Fateman. file directory \verb|http://www.cs.berkeley.edu/~fateman/generic|
%(http://www.sciencedirect.com/science/article/B6TYN-46C8RF8-1/2/cb147638a60fd464363ae5db7f337ede)
%Author Keywords: Subdivision; Interval analysis; Range analysis;
%Algebraic curves
\bibitem {kahan}
W. Kahan,
``How Futile are Mindless Assessments of Roundoff
in Floating-Point Computation?''
\verb|http://http.cs.berkeley.edu/~wkahan/Mindless.pdf|
\bibitem {berz} Hongkun Liang and Mark A. Stadtherr,
``Computation and Application of Taylor Polynomials
with Interval Remainder Bounds,''
{\tt http://citeseer.ist.psu.edu/cache/papers/cs/8022/ }
\bibitem{markst} http://www.nd.edu/~markst/la2000/slides272f.pdf
\bibitem {mpfi}MPFI group (Revol, Rouillier) INRIA
http://www.cs.utep.edu/interval-comp/interval.02/revo.pdf
\bibitem{nedialko}
Nedialko S. Nedialkov, Vladik Kreinovich, Scott A. Starks,
``Interval Arithmetic, Affine Arithmetic, Taylor Series Methods: Why,
What Next?'' (2003) \\
{\tt http://citeseer.ist.psu.edu/nedialkov03interval.html}
\bibitem {walster}G. W. Walster, ``Sun Forte Fortran,''
http://developers.sun.com/prodtech/cc/products/archive/whitepapers/tech-interval-final.pdf
\end{thebibliography}
\end{document}
;;;;;;;;;;;;;;;;; edited to here
\section{where is the payoff?}
\begin{verbatim}
{do I believe this? consider
$x(x+3)$ vs $x^2+3x$
Let x = [-2,1]
x(x+3) is
[-2,1]*[1,4] = [-8,4]
$x^2+3x$ is
[0,4]+[-6,3] = [-6,7], better, tighter..}
let x= [-2,1]
consider $x*(x+2)$ vs $x^2+2*x$.
[-2,1]*[0,3] = [-6,3] [0,4] + [-4,2] = [-4,6]
A tight bound would be [-1,3]. On [-2,1],
$f(w):=w^2+2*w$ has a min at w=-1, where f(-1)=-1
and a max at 1, where f(1)=3. Our simple routine for {\em NOT}
using Horner's rule
will not work for (x+A)*(x+B). We need to make
more use of symbolic info: the polynomial roots.
\section {How far to extend, e.g. beyond univariate polynomials?}
\begin{verbatim}
Ninterval^NInterval -> no change, or change to Sinterval.
We suspect In the following, k is the next unused interval index.
exp(Ninterval(a,b,#n,f)) -> Ninterval(exp(a){round down},exp(b){round up},#k,1)
log(Ninterval(a,b,#n,f)) -> Ninterval(log(a) etc)
sin(Ninterval(a,b)#n,f)) -> find if rel max min occur in [a,b], check ends
Ninterval(low (round down if not -1), high (round up if not +1) #k, 1)
cos (etc)
must fill out the matrix of programs for all binary ops and unary ops
Ninterval OP other
other OP Ninterval
* - * / expt
=, not=, >, >=, <=, <, included, overlap, disjoint.
\end{document}
.....................edited to here............
as bounded but not definite, or ``no number a`t all'' , or
``no interval at all'' as would be the
case for the empty intersection of two disjoint real
intervals\footnote
{One way of encoding such an interval might be [oo,-oo] which certainly
has nothing in its interior.}.
In some CAS designs, in particular Maple and Mathematica, the
designers have yielded to the temptation to have intervals
serve an additonal duty as representations for subsets of the real line.
For example, they compute $\lim_{x->\infty}\sin(x) = [-1,1]$. This is
a distinct interpretation of an interval and is known in the interval
arithmetic literature as a {\em cset}. By comparison
we remind the reader that ordinarily an interval is a
representation of a {\em particular scalar value somewhere on a real-line
segment}. We can't tell {\em which} point. The set interpretation of [-1,1]
as this limit is different, being {\em not} any such point, but
something like ``there is no particular limit but the function is
bounded
strictly between -1 and 1''. The usual definition of a limit requires
that as $x->\infty$ the distance between $r=\sin(x)$, which we know to
be $-1\le r \le 1$ and this computed limit, the interval $I=[-1,1]$ approaches
0, when in fact interval arithmetic tells us the distance is not zero,
but $[r-1,r+1]$. Since the same limit program (in Mathematica 10)
returns [-2,2] for the limit of $2 \cos x \sin x$, apparently (and
unsurprisingly) the boundaries need not be tight. In fact, $2 \cos x
\sin x= \sin (2x)$ and so tighter boundaries would be [-1,1]. Once
one admits loose boundaries, it seems it would not be ``incorrect'' to
return $[-\infty, \infty]$ for {\em any computation of a real
limit}. This is unlikely to receive full credit on a calculus exam.
% added 9/10/09 RJF , 9/30/09 RJF , 9/13/11 RJF
There are further ramifications: computing $\sin^2\infty$ as [0,1] in
Mathematica suggests that in some obscured form we might start with
something equivalent to $\sin^2x+\cos^2x$ but a failure to
sufficiently simplify followed by an evaluation at $x=\infty$
reduces this to $[0,2]$. This may be acceptable in some circumstances,
(like if you find it acceptable for $\lim_{x \rightarrow \infty} 1 ~=~ [0,2]$.)
but it points to a need for the user or the program controlling the
computational flow to check for intervals or unconventional
operands pervasively, {\em and not only whenever limit calculations are used.}
This may slow down ordinary computations on ordinary objects, for
which such shortcuts as ``$0$ times anything is $0$.'' almost always prevail.
Indeed, using such ``generically correct but sometimes wrong''
rules seems unavoidable, the only proposed remedy has been
to carry along with such a computation some
set of provisos like ``$0$ unless $x$ is a $\bot$'', or a rigid
context-setting environment in which any operands or operations
which cannot be validated cause an error exception in the flow of
a procedure.
%Limit[Sin[x+Pi/4],x->Infinity] gives [-sqrt(2),sqrt(2)] in ver 9,
%[-1,1] in version 10... note ...
%$$\sin \left(x+{{\pi}\over{4}}\right)={{\sqrt{2}\,\left(\sin x+\cos x
% \right)}\over{2}}$$
A more pervasive issue is the alteration of the limit program from an
equivalence transformation (ET) to an approximation program (AP). By
an ET we mean a program that changes one
form\footnote {A more careful
formulation would worry about domains. Here we are, with some
trepidation, asserting the domain is ``whatever the CAS knows
about''.}
to another without changing its value ({\em perhaps}) with
some singular exceptions. Ordinarily the idea behind an ET is to
provide an equivalent but, by some measure, simplified expression.
Thus an ET simplification program might reduce ${{2}\over{\sqrt{2}}}$
to $\sqrt{2}$ or $x^2/x$ to $x$. In the latter case we have found a
form which is the same\footnote{except when $x=0$ (arguably).}. We
expect all routine simplification transformations to be ETs, although
we may be occasionally disappointed when they are not. This may
happen routinely when an expression depends on some implicit choice,
say of branch cut -- the simplest cases probably the consequence
of dropping the sign in an expression like $\pm\sqrt{1}$ or trying
to choose the ``right'' branch of $\log(x)$ knowing nothing about
the use being made of it. In particular, the value of $x$ does
not determine the right branch..
On the other
hand, an AP is a program that given an expression, can produce
\begin{itemize}
\item (a) an approximate numeric answer. The request
to obtain such a numeric result is often signalled implicitly by
inserting a floating-point number as in
$\sin(3.0)$ which may result in the approximation
0.14112000805987 or
\item (b) an explicit request for a numerical
valuation, or
\item (c) a symbolic approximation such
as a request for a Taylor series approximation:
$$\sqrt{x+1}=1+{{x}\over{2}}-{{x^2}\over{8}}+{{x^3}\over{16}}-{{5\,x^
4}\over{128}}+\cdots $$
\ene{itemize}
While no one legislates such matters, we expect that
ETs follow certain rules, like ET(a)-ET(b) is equivalent mathematically
to some ET(a-b), assuming appropriate domains.
Ordinarily a limit calculation would be an ET. This assumption is
violated if it returns $\bot$ or an interval or $\pm \infty$.
Arguably, then, if we
want limit() to be an ET, it should return itself, that
is, ``unevaluated'' when it would otherwise be $\bot$ or an
interval. This would allow certain
calculations to proceed by (for example) allowing limit expressions
to contain their original information in such a way that
they might be combined with (or cancelled
by) other related limits which might otherwise also have infinite values.
As an aside, it is possible to group other commands in a CAS into
two more classes:
\begin{itemize}
\item Valuation transformations (VT) in which an expression is
evaluated by substitution of values for indeterminates, e.g. $x=1$.
Numerical approximation could be viewed as a VT where the values are
approximate numbers, but some valuations are exact or even symbolic as
$x=\cos(\theta)$.
\item Derived transformations (DT) in which some new expression
(non-equivalent) is derived, perhaps by differentiation,
integration, or simply by applying arithmetic operations and other
combinations to expressions.
\end{itemize}
(There are also input and output commands for reading, writing, plotting, etc.)
%***** 10/1/2009 RJF end
\medskip
\section{A closer look at intervals}
Let us revisit that assertion that $[-1,1]-[-1,1]$ might be zero, or
not. This question is referred to as the dependency issue, in interval arithmetic
(or ``reliable computing'') circles, in which the inability to
determine the ``source'' of an interval, i.e. whether two intervals
are correlated, leads to wider and therefore less useful results.
In our example, if the expression means the set $\{x-x ~|~ x \in
[-1,1] \}$ the answer is 0. If we mean \{$x-y ~|~ x \in [-1,1]$ and
$y \in [-1,1]$\}, then the answer is the interval [-2,2].
Is there a neat way of distinguishing the two cases? Here's a possibility,
using the CAS Macsyma, where these might be notated
\begin{verbatim}
block([x:interval(-1,1)], x-x)
\end{verbatim}
and
\begin{verbatim}
block([x:interval(-1,1),
y:interval(-1,1)], x-y)
\end{verbatim}
This locution is an attempt to force a kind of interchange of
simplification and evaluation, where in the first case the expression
{\tt x-x} would be simplified to zero and in the second case the
expression x-y would not be simplified, but just evaluated using
interval arithmetic rules\footnote{Without changes in Macsyma, which
currently does not have rules for arithmetic about intervals, it
doesn't work. They both return 0.}.
Another example which illustrates the possibility of getting tighter
bounds with additional work, is the expression $\sin(x) \cos(x)$
which, on the face of it, for all real $x$, seems to be in the
interval [-1,1]. In fact, the expression is in [-1/2,1/2], a result
which becomes clear by recalling that $\sin(x) \cos(x) = \sin(2x)/2$.
A non-constructive definition of interval arithmetic can be phrased in
terms of constraints on a function. Elucidating the situation of
dependency requires more than one variable, so let us start with two.
Then given a real-valued function of two real-valued arguments to
evaluate: $f(x,y)$ for $x,~y$ interval-valued between $[x_0,x_1]$ and
$[y_0,y_1]$ respectively, we return an interval $[z_0,z_1]$, where the
real-valued $z_0$ is the minimum of $f(x,y)$ with the constraints $x_0
\le x \le x_1$ and $y_0 \le y \le y_1$, and the corresponding maximum is
$z_1$. While the interval $[z_0,z_1]$ then is the ``best'' answer, it
might not be possible to find it conveniently either because we have
lost information during the calculation, or the minimization
computation is computationally too costly. We tend to accept the
computation of some other, looser, version: $[w_0,w_1] \subseteq
[z_0,z_1]$. That is $w_0 \le z_0$ and $w_1 \ge z_1$. We may also need
to extend the interval notation to allow for situations such as division by
zero.
%%\begin{verbatim}
%%\end{verbatim}
Much effort can sometimes be expended in improving ``tightness'' of
the result, and this may be useful in the long run, since any slack in
the result can be magnified by subsequent operations.
In our assumptions about a CAS, we mean by ``interval()'' a program
which produces a new object, not shared.\footnote{Maple has a feature
making all compound objects which simplify to the same structure
isomorphic. Implementationally they are pointers to one copy in
memory. This is helpful in some processing, but not this one. There
is sometimes a difference between holding two chocolate-chip cookies
behing your back,
one in each hand, and holding just one cookie but with both hands.}
\section{An approach}
In this section we discuss a collection of approaches for computing
with {\em real intervals}. We claim no novelty in these
matters\footnote{Indeed, any time we thought we had a new idea in
interval arithmetic, we found we were simply rediscovering something
that had been previously identified as interesting (and often analyzed
in great detail).}. Our perspective here is to briefly introduce
ideas to a reader who may not have seen these before, or has not seen
the alternative possibilities to `naive' intervals. Our objective is
to provide enough background so the reader can seek out further
information in the references, but also so that we can reasonably, in
the following section, point out the relationship with pseudo-values
that occur in CAS.
\subsection{Intervals are Three-tuples}
To be upfront about this approach: it is not one that we recommend
in the context of a computer algebra system (CAS). In a very real sense
this
approach effectly inserts a mini-CAS into the interval library.
We create a numeric interval of the real line by calling a function,
{\tt Ninterval}. For example, for the interval between -1 and 2, we
evaluate a form like this $v=${\tt Ninterval(-1,2)}.
%(setf v (ri -1 2))
You might expect that this function would construct something
reflecting the usual representation of an interval as a pair,
e.g. [-1,2].
A slight variation is used: we produce
an internal data structure with {\em three} slots. The
% (describe *)
%[-1,2] is a structure of type ri. It has these slots:
% lo -1
% hi 2
% exp #(#:g478 0 1) ;; note, g478 is a made-up index.
first two are the visible upper and lower limits, sometimes denoted
${\underline x}$ and ${\overline x}$, where in our example, ${\underline
x}=-1$ and ${\overline x}= 2$. These can be exact integer or rational
numbers, machine floating-point numbers, software ``bigfloats'', or
(though we do not recommend this) literal parameters or expressions (say $a$
or $\exp(b-1)$). The presence of literal parameters should result in
almost all operations on such intervals resulting in a non-evaluating
``quoting'' of the operation.
Our first proposal for the third slot is to manufacture a new symbolic
variable name\footnote{a counter is incremented after each newly
constructed name. There are complications if one tries to merge
saved results from two otherwise unrelated computational sessions:
coincidental counters must be relabeled.}.
Here we denote that
variable by \#n for some integer n. The third slot is then an
algebraic expression involving that variable, initially just \#n. The
exact representation used is specified in the implementation, and
might be a list of polynomial coefficients, a ``tree'' or some other
form that might be convenient in a CAS.
Our example $v$ might be described this way: {\tt[-1,2,
\#5]}. However, since we (humans) ordinarily do not need to see the
extra slot, programs ordinarily will not display it in the computer
output. For purposes of exposition in this section of this paper, we
will show them.
If we compute $u= v \times v$, we create the slots in $u$ as
{\tt[0,4,\#$5^2$].} Note that if we had two independent
intervals that happened to have the same bounds, and multiplied them,
we would compute {\tt[-1,2,\#5] * [-1,2,\#6]} as {\tt [-2,4,\#7]}.
This is the interval [-2,4], which is wider than [0,4]. Furthermore
we have generated a new variable, {\tt \#7} which means we
have lost some information about the origin and possible restrictions
on the values of
this interval, if, say, we were to further combine it with {\tt \#5}
or {\tt \#6}.
Let us continue our example.
If we compute $u \times v$, we can notice a coincidence in indexes in
{\tt [-1,2,\#5] * [0,4,\#$5^2$]}. By realizing that we are
computing \#$5^3$, the resulting interval is [-1,8], instead of the result of straight
multiplication of intervals, which would give [-4,8].
{Under certain conditions we can maintain this stream of computation,
and more generally recompute the expression in the third slot. To do
this we must not only represent the polynomial, but remember the
original value. In this example we do this by keeping the original
{\tt \#5} around as long as there is a live ``accessible'' value that uses
it. There is a convenient mechanism that is built-in to Lisp for this
purpose, a ``weak hash table''.}
A bit of introspection as to what might be appropriate for the third
slot
suggests that further pushing this idea results in building a CAS to
manipulate it. It is probably better to do the manipulation needed
with indeterminates and use the more-or-less full power of an existing
CAS if the environment provides it. What might this look like?
%11/13/2015
\subsection{Intervals are Pairs}
Here's a simpler approach which takes fuller advantage of the CAS, an
approach in which the desired optimization of interval arithmetic
is obtained by deferring the ``interval evaluation'' to a stage
after the symbolic simplification or rearrangement of the expression
to be calculated. To follow
the same example as above,
If we compute $u= v \times v$, we defer any interval computation
until we simplify the expression to $u=v^2$. If this is as far as we
need to go, we can evaluate to an interval by computing
{\tt power([-1,2],2)}.
If this is not the final stage, but we next compute $u \times v$, we
learn that (symbolically) we have created $v^3$. This can be
evaluated as {\tt power([-1,2],3)}. This approach ``externalizes''
the dependencies rather than trying to keep them internal to the
third slot.
How to deal with {\tt v=[-1,2], v=v*v}? It is simpler to
deal with a computation which is essential ``functional'' in
nature, but assignments can be rephrased as (essentially)
the generation of a sequence of values for {\tt v}
where we separate the current symbol expression from the previous
one(s) and ultimately link these to the value of {\tt v} from the
original statement.
%11/13/2015
%Another question or a few: What would programs look like?
For example, to evaluate an expression
where {\tt ri(a,b)} constructs a real interval, {\tt widen\_ri(c)}
constructs the smallest interval containing c:
\begin{verbatim}
intervalbind(
[v: ri(-1,2),
w: widen_ri(3),
u: ri(sin(3.0d0)),
v2: ri(-1,2),
x: 1/2]
[ v^2-1, v*v2-1, (v+1)*(v-1) ,
w*w, w^2, v+x, u-sin(3.0),
v*interval(-1,2)])
\end{verbatim}
What would an iterative interval program look like?
We cite Nathalie Revol
http://perso.ens-lyon.fr/nathalie.revol/publis/RR05.pdf
A simple version of an Interval Newton method might add to
the usual iteration, a check that the interval iterate
contracts around a zero, and if it does not by a Newton step, it
reverts to a bisection.
A well-known optimization we can make symbolically
applies to any quadratic polynomial. We can
rewrite $p=3x^2+2x+1$ as $q=3(x+1/3)^2+2/9)$. This is a ``single use
expression'' (SUE) \cite{walster} referring just once to $x$, and
provides tighter bounds as $p([-1.1,-1.0])$ is $[1.8,2.63]$ but
$q([-1.1,1.0])$ returns the tighter (in fact optimal!) $[2.0, 2.43].$
The order of computation here is based on ``completing the square.''
Unfortunately this particular trick is specific to quadratics, and
generalizations are not as effective. We can reduce the number of
appearances of $x$ by rearranging terms, e.g. for a cubic we can
arrange for only two appearances of $x$. First divide through by the
leading coefficient so that we have $x^3+r x^2+s x+t$. Then by
setting $(a+x)^3+b x^2+c= x^3+r x^2+s x+t$, we can find values for
$a,~b,~c$. Let $v=\sqrt{s}$, then $\left\{b = r-\sqrt{3} v,c =
t-\frac{v^3}{3 \sqrt{3}},a = \frac{v}{\sqrt{3}}\right\}$
We have experimented with some of the theoretically faster ways of
evaluating a polynomial via rational preconditioning calculations with
the coefficients, but these seem to fail to produce an improvement. A
much more complicated procedure that (at compile-time)
essentially divides up a polynomial $p$ with scalar coefficients into
monotonic segments $\{p_i\]$. Then the upper and lower bounds for each $p_i$
can be computed by scalar computations at the endpoints of
the interval. The results of $p\i({\underline x})$
and $p\i({\overrline x})$ will generally result in intervals, again.
These individual bounds on pieces have to be sorted for a global
result.
Tighter results might be obtained by evaluating $p_i$ in extra precision.
Of more interest are multivariate
expressions, when nice tricks are harder to come by. \cite{demmel}
%http://www.cs.berkeley.edu/~demmel/demmel_talk_Banff2005.pdf
A tangent regarding implementation in Maxima. This CAS, implemented
in Lisp, provides an internal list structure for storing properties.
For example, $x^2-1$ is {\tt ((mplus simp) -1((mexpt simp) x 2))}
where the {\tt simp} indicators are used to mark the subexpressions
as already simplified. This complicates the structure, which
conceptually might be written as {\tt(mplus -1 (mexpt x 2))}, but has
an important use in avoiding redundant simplification computations.
This encoding for attaching properties to expressions has been
in use since the initial design, circa 1966.
Rather few uses have emerged\footnote{ In retrospect, if we knew
so few properties would be developed over the years, other schemes may come
to mind. For instance, if there is one property ``simplified'' then
all such expressions can be implicitly marked by being stored in
particular machine addresses above some boundary. If there is a
property
``evaluated more recently than expression X'' then expressions can be
stored in tagged pages of memory. Both these techniques are used in
some other existing CAS systems.}:
here's an illustration of two uses.
If we factor that expression to yield $(x-1)(x+1)$, it is
marked up this way: {\tt ((mtimes simp factored))((mplus simp
irreducible) -1 x)((mplus simp irreducible) 1 x))}.
In the interval encoding, it is tempting to put the third, index, argument in as a
``property''
for Maxima, making interval(-1,1) look internally like
{\tt ((interval simp g478) -1 1)} or perhaps
{\tt ((interval simp (index g478)) -1 1)}.
After some experimentation it becomes clear that this would not have
all the right properties because Maxima relies on an assumption that
is, in this case, violated. The assumption is that if two
general form expressions X, Y are identical
without looking at the properties, then X+Y is 2*X or perhaps 2*Y.
That is, the properties are to be thought of as
not specific to that object in memory,
but are presumed true for any isomorphic object, where isomorphism
compares ``heads'' like {\tt mplus} or {\tt interval} but not the
properties.
Consequently it is not a workable representation in Maxima to
use the property list. We proceed to use a third argument as
an index, or (see the next section) some other information.
%*** changed 11/13/2015
\subsection{Variations}
In the absence of a CAS framework, the third slot becomes
more appealing. It is tempting to write a mini-CAS to deal
with the symbolic information about
the interval that could be internalized in that third slot,
What should we allow as acceptable expressions for that third slot?
We initially considered just expressions of the form $r x$ for a
scalar number $r$ or $rx^n$ for some integer $n$ but then realized
that it would be simple to use existing subroutines for symbolic
mathematics to carry around a more elaborate form. The next
elaboration was a polynomial in (any) {\em one variable}, in some
encoded canonical form. {Given such a form, we can feed it into a
program for further analysis: for evaluating a polynomial $P$ in one
variable $x$ at an interval point, the tightest bounds can be
computed by considerably more elaborate computation than just going
through the adds and multiplies: locating suitable intervals around
the zeros of $P'(x)$, then evaluating $P$ at those places to find
relative extrema. Combined with $P$ evaluated at the endpoints (absolute
min/max), we can find sharp bounds subject only to the representation
of the endpoints. Such a computation {\em could} be worthwhile as
a trade-off: the advantage of finding tighter bounds versus the
considerably increased computational expense.}
Any third-slot expression that turns out to be more complicated than we wish to
deal with can simply be given a fresh name, and from that point forward it is
assumed to be independent of all previous variables. Thus our current
program would never notice the correlation possibility that is
inherent in the identity $\sin(x) \times \cos(x) = \sin(2 x)/2$ since
each of the expressions $\sin()$ and $\cos()$ would be given different
(presumed to be independent) names. This would lead to more conservative
bounds, but not incorrect ones.
\medskip
Where do we draw the line, though: what is too complicated? Can we in
fact use something more elaborate than polynomials? Indeed, there are
a host of variations on interval arithmetic: we found that our invention of a
``polynomial'' slot is subsumed under any of a number of previously
documented innovations, many of which are reviewed by Arnold Neumaier,
see \verb|http://www.mat.univie.ac.at/~neum/papers.html#taylor|. This
paper reviews so-called Taylor forms---an extension which revolves
around carrying a truncated Taylor series with remainder, evaluated at
the midpoint of the interval. In certain calculations this additional
information can be effective in reducing the width of the intervals
being carried. To produce this model, we could use the extra slot
effectly as an expression of a Taylor series, with remainder, for the
number being represented. We could start our trigonometric
calculation by considering that we are evaluating $f(x)=\sin(x)$,
where $x=[a,b]$. We then expand sin in a Taylor series around the
midpoint, $x=x_0=(a+b)/2$, and use an expression for the remainder
(error). Most CAS provide a suite of operations on Taylor series, and
so can implement this ``off the shelf''. The operations are similar to
those on polynomials with the major difference that the Taylor method
sets a global order, (maximum degree) and high-order terms are
truncated. The error term is calculated separately as another interval
calculation.
%[[[Just refer to Berz paper in Reliable computing...
It is important to observe that the evaluation of the remainder term,
which can itself be elaborate, may have to be done carefully too, or
some blow-up is re-inserted. Programs for evaluating derivatives can
be generated by ``Automatic Differentiation'' programs, but these must
also be evaluated in an interval fashion. When the error term looks
sufficiently well-behaved, ``naive'' interval arithmetic may
do. Illustrations of the payoff from maintaining tighter bounds seem
tied to a procedure which can terminate sooner because of the bounds.
A typical case may be a Newton iteration combined with interval
bisection (or some subdivision, since an interval with an infinite
endpoint cannot really be bisected) which may converge with {\em far
fewer iterations} and thereby pay for the more expensive
operations.
%]]]
Other possibilities include allowing the extra slot to contain a
linear function of any number of interval values (more about this
later). In this case a linear-programming method could be used to find
min and max. If the extra slot were allowed to be a multivariate
polynomial, a method based on cylindrical algebraic decomposition
might be suitable.
An idea called ``Affine Arithmetic'' for intervals has also been
pursued in the recent literature on reliable computation. In this
model, a quantity $x$ is represented as a first-degree (``affine'')
polynomial in the values it depends upon. In particular, $$x=x_0+
x_1e_1+x_2e_2+ \cdots + x_k e_k$$ where the $\{x_i\}$ are given ``real
numbers'' and the $\{e_i\}$ are dummy variables, assumed uncorrelated
for calculation purposes, but each of whose values is known to be in
$[-1,1]$. When possible the representation is maintained through
ordinary arithmetic, and when it is not, the operational routine must
find an alternative acceptable approximation by introducing another
variable $e_{n}$ and positing it as uncorrelated. Affine arithmetic
(AA) seems especially appealing when geometric objects are defined in
which the coordinates share a set of dummy variables. This situation
produces centrally-symmetric convex polygons or ``zonotopes''. In two
dimensions these can be visualized as ``lozenge'' shapes that are
smaller than the enclosing rectangular intervals that would appear
from naive intervals. There are benchmarks illustrating particularly
advantageous results in approximate computations for graphics: they
can provide closer adherence to lines or planes that are being
approximated. Not all computations ultimately benefit from AA:
it may be beneficial to switch between methods to
increase overall efficiency: AA may be better for small
intervals. An extensive comparison of interval variations is given in
the paper by Martin {\em et al}.\cite{martin}.
%(a +b e1) +( c +d e1) --> (a+c) + (b+d)e1.
%(a +b e1) *( c +d e1) --> (ac) + (ad+bc)e1 +bd*e1^2
% A*e1^2 -> A/2 +A/2*e3
See for example,\\
\verb|http://www.icc.unicamp.br/~stolfi/EXPORT/projects/affine-arith/Welcome.html|
or\\
\verb|http://www.tecgraf.puc-rio.br/~lhf/ftp/doc/oral/aa.pdf|
%http://www.ic.unicamp.br/~stolfi/EXPORT/projects/affine-arith/Welcome.html
Implementing AA on top of our structure is straightforward, and we
have done so up to addition and (non-trivial) multiplication, for
arithmetic combining scalars and intervals.
These ideas are aided by such tricks as the previously
mentioned ``single use expression'' (SUE) from some set of variables\cite{walster}. This would transform an expression like $ax+bx$ where $a$ and $b$
are scalars, but $x$ is an interval, into $(a+b)x$, which may be a
tighter expression.
Similarly, for intervals $x$ and $y$, $A= x/(x+y)$ can be rewritten as
the SUE $B=1/ (1+y/x)$, although the acceptable domains for $x$ and
$y$ have now been changed. (Look for division by interval expressions
that might include zero!) Walster suggests obtaining the answer by
computing $C = 1- (1/(1+x/y))$ and then returning their intersection,
$B \cap C$.
One approach for SUE is for us to establish a set of heuristics for
driving an algebraic expression into a more-nearly SUE form. A good
algorithm to find an optimal ``sharp'' solution (that is, more
efficient than exhaustive enumeration and evaluation) would be
welcome, but tradeoffs such as whether a SUE version with respect to
$x$ is better than a SUE with respect to $y$, or if some other
combination is yet better, cannot be computed at compile-time, since
in general the optimimum depends not only on the expression (which is
presumably available at compile-time), but on the particular {\em
run-time} interval values associated with $x$ and $y$.
{A possible, but somewhat unappealing, approach can be constructed
around multiple arrangements of the same expression. Given any
collection of ``mathematically equivalent'' expressions transformed at
compile-time into various interval expressions, each version could be
numerically evaluated, and the tightest solution emerges by computing
their intersection.} Such evaluations must be done while aware of the
possibility of the introduction of extraneous singularities; if any of
the expressions is non-singular at a point in the interval, then it
presumably reflects a valid mathematical result. A good expression
might be a Taylor series expansion at the midpoint of the
interval. This cannot be computed without knowing the interval
endpoints. An interesting approach shows that an endpoint of $f(x)$ can be
shown to be sharp if it depends only on one endpoint value of the input
interval $x$ \cite{hansen}.
We could attempt to use some existing (essentially black-box and
possibly heuristic) programs like Mathematica's {\tt Maximize} and
{\tt Minimize}. Yet, even when they nominally work according to their
own specifications, the answers may not be satisfactory for interval
computations. Reliable computing requires strict inclusion: these
imply stronger requirements for directed rounding of endpoints than
typical numerical min/max programs.
In Nedialko {\em et al} \cite{nedialko}, various methods are
considered. A theoretical comparison can be summarized by an
optimality result: representation by polynomials is, in their
presented context, optimal.
\subsection{Other approaches to improved interval computation}
Other techniques, such as (repeated) bisection of intervals can be
used to attempt to narrow the result width. Our suite of programs
includes this conceptually simple technique: the result for $P(x)$
would be the smallest interval containing all the ranges of $P(x_1),
..., P(x_n)$ where $x$ has been subdivided. That is, $x$ is the union
of the $x_i$; some criteria determine which of the subdivisions to
divide again in an attempt to narrow the interval. (In particular the subdivisions
containing the current maximum and minimum should be re-divided.)
This kind of interval arithmetic is more encapsulated and therefore
``hands-off'' compared to the usual numeric version. It takes the
form of a functional programming model: new values may be created out
of components; but we do not allow user-level operations like ``change
the lower limit of variable $x=[{\underline x},~{\overline x}]$ to be
the number 43.'' Instead, we allow the creation of a new object which
can be re-assigned to the same name: $x=[43,~{\overline x}]$. Old
values, if they are no longer accessible, will be collected, and their
memory locations re-used.
\section {Pseudo-values: Infinities and intervals}
At this point it may not seem that there is any direct connection
between intervals and extended numbers such as infinity. In fact,
the resemblance has more to do with the techniques to resolve
difficulties than the mathematical concepts.
Note that without intervals, the production and therefore occurrence
of infinities is limited to singularities (of which division by zero
is the most obvious, but not the sole example). With intervals, the
exceptional treatment of zero-divisors must be extended in some way to {\em
intervals that contain zero}. Once having produced such an object,
can one make it go away? For example, can we be assured that $1/(1/x)=x$
true, even if $x$ is an interval containing zero?
\subsection{Infinity as a symbol}
One of the problems of a CAS design is that there is a tradeoff
between convenience for human users and for computation. As an
example, computationally, we might construct a system in which an
infinite value appears only as a result of certain operations, and in
each case we can find some more informative alternative. Unfortunately,
if its origin is simply ``a human typed it into the system'' we do not
have an alternative to ``just compute with this.''
What can we propose? We can keep $\infty - \infty$ from becoming zero
if, each time we generate (a new) $\infty$, we give it a unique
subscript or index, in a manner similar to the indexed variables or third slots used
with the intervals. That is, $\infty_1$, $\infty_2$ ,..., $\infty_n$.
In Macsyma/Maxima this is rather simple to start:
{\tt (inf\_counter:0,
infinity():=(?incf(inf\_counter),Infty[inf\_counter]), nofix(infinity))}
After this, any mention of {\tt infinity} as in {\tt infinity+3} will generate
an expression like {\tt Infty[4]+3}. In \TeX\ this could be displayed
as
$\infty_4+3$.
This is not enough: A search through the Maxima code for all producers
and consumers of infinities may reveal locations that must be changed
to conform to this model. These occur principally in the {\tt limit},
{\tt sum}, and definite integration facilities. Also saving such
values off-line and reading them back must re-create the proper
values, shared only when appropriate, and generally distinguished from
values created from other processes and possibly also saved.
Maxima distinguishes between
real infinity({\tt inf}) and complex infinity ({\tt
infinity}). Mathematica 7 and later correspondingly uses {\tt
Infinity}, internally {\tt DirectedInfinity[1]} and {\tt
ComplexInfinity} actually {\tt DirectedInfinity[]}. Other directions
are possible, e.g. {\tt DirectedInfinity[Complex[0,1]]}, where the
direction is the vector from the origin to
a designated point on the unit circle.
$(1+i)\infty$ is effectively {\tt DirectedInfinity}($(1+i)/\sqrt{2}$) after
normalization to unit length.
A set of indexed infinities in
each direction is a straightforward elaboration of the ideas
presented here.
We can also link the directedness to the representation, and so we
will subsequently consider DirectedInfinity(direction, index).
We need a new simplifier that is able to compute a kind of pessimistic
simplified form. For example, given the expression $\infty_2 -
\infty_3$, {\em assuming that there are no occurrences of} $\infty_2$
or $\infty_3$ {\em anywhere accessible in the CAS}, can be simplified
by generating an undefined: {\tt und}. Such symbols are also
indexed. If there are other occurrences of those $\infty_k$ symbols
elsewhere, it is prudent to keep the expression uncombined in the hope
that at some later point one of those symbols will interact and
cancel. We typeset the undefined quantity with the ``bottom'' symbol,
for example, $\bot_4$. This gives us an opportunity to consider a
case where we assign {\verb|x := limit(1/t^2,t->0,plus)|} and thus
x=$\infty_5$ , a specific value. Then {\tt x-x} is $\infty_5-\infty_5$
, and indeed this can be simplified to $0$ because this is the
difference of two of ``the same $\infty$''. In a system with directed
infinities $-\infty_5$ would have to be correlated with
$-\infty[-1]_5$ for such cancellation. Mathematica's current scheme
presumably has only one infinity in each
direction, and so rules combining them must either be very
conservative or sometimes wrong. This design suggests that for each
directed infinity generated, there should be a subscript relating to
the magnitude of the infinity as well as a separate direction.
This simplifier must deal with any expression with one or more
infinite or undefined objects, possibly other forms such as
``indefinite'' (but bounded) and generally should come up with the
most informative and compact result possible. For example, {\em in
the absence of any other occurrences of the same-indexed infinities,}
$\infty_2+\infty_3$ could be replaced with $\infty_4$ as a simpler
result\footnote{ To reiterate: it would be a mistake to make this
replacement if there were another expression algebraically dependent on
either of the two component $\infty$ objects, since the replacement would
lose dependency information.}. In similar circumstances (no
dependencies) there is no point in retaining the details of
$(\bot_4+1)/\bot_7$. This should be replaced by (say) $\bot_8$. An
expression involving such objects should be simplified as much as
possible, using standard rules for combination.
This is not necessarily trivial, and for rational functions
may be solved using a resultant technique devised by Neil Soiffer in
his MS project at Berkeley (1980). That work calculated the minimal
number of ``arbitrary constants'' needed to provide the same degrees of
freedom as a given expression with perhaps many such constants.
%%% new section
\section{Arithmetic with Pseudo-values }
The failure of elements in the classes of $\infty$ and $\bot$ to
satisfy some of the ordinary axioms of real or complex numbers has two
possible consequences: either this invalidates certain standard rules
of {\bf operation} in a CAS or limits the mathematical {\bf domain} of the CAS to
exclude pseudo-values as a full-fledged participant in computation.
Since the ambition of at least some CAS is to encompass all of
mathematics, such systems must either entirely march $\bot$ out of
mathematics, restrict such results as ``output only'' say as something
that can be printed, but not suitable as input to any further
operation, or figure out how to ``do something''.
Examples of potential failure
%*****
zero identity $0*x=0$. What if x=inf or und?
multiplicative inverse $x* (x^(-1))=1$. What if x=inf or und?
cancellation if $a*x=b*x$ means $a=b$, does $3*inf=4*inf$? and does that mean 3=4?
$x+y>x$ if $y>0$. what if x=inf?
is inf=inf? und=und? what about NaNs?
How can we address this? One way, which seems impractical, is to
change
the simplifer so that 0*anything is simplified to
``0 unless anything is in {und or inf}''.
This kind of information must be carried along throughout
the system, and can be very clumsy. If one can prove that
anything is (say) explicitly a finite constant like 43,
then the ``unless'' clause can be dropped. Otherwise it has
to be carried along, and combined logically with other ``unless''
clauses.
What to do then?
One choice. If you have a zero in the algebraic real field within which we
are usually computing, 0* anything is 0. End of discussion (and there
is at least some agreement on this) even if the anything is part of
the extension, i.e. oo or -oo.
Kolmogorov, N. A. "Infinity." Encyclopaedia of Mathematics: An Updated
and Annotated Translation of the Soviet "Mathematical Encyclopaedia,"
2nd ed., Vol. 3. (Managing Ed. M. Hazewinkel). Dordrecht, Netherlands:
Reidel, 1995. {page 193} cited in
http://mathworld.wolfram.com/AffinelyExtendedRealNumbers.html
Here's a modification: If that is not zero, really, but
z=1/oo\_i, then oo\_j*z is indeterminate unless i=j, and then
it is 1.
Maybe make it illegal to do certain assignments, namely
giving a name to a quantity that
is und or inf.
For example {\tt z : limit(1/x^2,x,0)} is an illegal {\bf assignment.}
(That is, inf is essentially ``He who must not be named'').
We hope to make it impossible to encounter a situation of 0*z in which z
subsequently is assigned infinity.
There is also the need to disallow substitution as subst(inf,z,0*z)
which would be too quick to simplify 0*z to 0.
This also eliminates the possibility of z-z being inf-inf, or
0*z being (through a later assignment to z) being 0*inf.
For brevity let us define
infmaker():=limit(1/x,x,0)
We do not allow z:infmaker()
but what about z:1/infmaker() ? Is this 0? It appears so. Then
even if we have uniquely indexed oo objects, 1/oo[i] would
be zero.
1/infmaker() * infmaker() is 0*oo[j] and we are left with the
same dilemma, without any assignment at all.
This computation o*oo should be $\bot$ not 0.
*****
oof.
How about this. We scan an entire expression and reserve judgment
on products (and sums?) and determine if any of them have components
which are inf or und (or +- inf, und). If so tread very carefully
.. essentially call a special evaluator... This still requires that
there be no assignments z:inf in the future..
Or -- for each simplification of mtimes or mplus, provide a
catch/throw
sort of like (catch 'thereisinf (simptimes ...) (recompute it...))
and within simptimes --- (dynamic return) if there is an inf or
und, do a (throw 'thereisinf)
Still remaining problem is that if there is ((mtimes) 0
)
then it seems that there is a prospect that the is not
examined at all!
This ``short-cut'' simplification either is present in Maxima now.
If it was implemented in the past, it has been disabled.
% proof
% kill(a,b,z)$
% tellsimp(a(3),(print(a3),0))$
% tellsimp(b(3),(print(b3),0))$
% tellsimp(z(3),(print(z3),0))$
% 0*a(3)*0*b(3)*z(3)*0;
%tellsimp(makeinf(1), throw(ReEvaluateMe)$
%catch (a*makeinf(1));
% returns ReEvaluateMe
%Write a program that simplifies mtimes (number, infContainer,
%symbols).
%If infContainer is (inf-inf) then und, if inf^(neg) then zero, if
%something else??
%****
%matchdeclare(any,any)$
%tellsimpafter(inc[any], throw (gottaInf));
%infsimp(r)::= block([res1: catch(ev(r))],
% if gottaInf=res1 then
% block([simp:off] ,?print(r),helpme) /*Here's where we fix up simplif.*/
% else res1);
See program in tex comments.
This section raises questions, but does not
answer all of them. The approach taken by a particular CAS may be
delimited by the activities it is supposed to support. Thus a base
for a theorem-proving system might be more ``pedantic'' than one
intended for (say) checking hand computations where it is known
beforehand that extraneous solutions will occur and can be discarded.
We start by observing the properties that are consistent and those that fail when we have $\bot$ and
$\infty$ without subscripts. (or in some cases, with subscripts).
Consider that at the base of the ordinary simplifier we assume that
the otherwise uninterpreted symbols $x$, $y$, $z$ etc., come from the properties
of integral domains, which include the postulates for commutative rings.
{\bf Commutative Rings:}
\begin{enumerate}
\item %1.
Closure:
$a,b$ in R imply $a+b$ and $a\times b$ are in $R$.
\item %2.
Uniqueness: $a=a'$ and $b=b'$ imply $a+b=a'+b'$ and $a\times b
=a'\times b'$
\item %3.
Commutation: $a+b=b+a;$ $a\times b=b\times a$.
\item %4.
Association: $a+(b+c)=(a+b)+c$; $a\times(b\times c)=(a\times b)\times c$.
\item %5.
Distribution: $a\times (b+c)=a\times b+a\times c$.
\item %6.
Additive Identity (Zero). there is an element, 0, such that $a+0=a$ for all $a$ in $R$.
\item %7.
Multiplicative Identity (One). there is an element, 1, such that $a\times 1=a$ for all $a$ in $R$.
\item %8.
Additive Inverse: $a+x=0$
has a solution $x$ in $R$ for all $a$ in
$R$. ({\em This would fail for $a= \bot$ or $a=\inf$ but is solved by
subscripting.})
%\end{enumerate}
{\bf Integral domain adds these}
%\begin{enumerate}
%9.
\item Cancellation: if $c$ is not zero and $c \times a=c \times b$,
then $a=b$. ({\em fails for $c=\bot$ or $c=\infty$}. We could
restate so that $c$ is defined and finite, but then we must place
this condition on every denominator.)
Other property usually assumed for the symbols include {\em ordering.}
%10.
\item Trichotomy: $a>0$, $a<0$ or $a=0$. ({\em fails for $\bot$})
%11.
\item Transitivity: $a****b$. ({\em fails for $\bot$ and $\infty$})
{\bf What about fields, e.g. reals?}
%15.
\item
A field is an integral domain in which each element $e$ except 0 has an
inverse, $e^{-1}$ such that $e\times e^{(-1})=1$. (This seems simple enough, but has lots of further
consequences.) ({\em fails for $\bot$ and $\infty$})
)
\end{enumerate}
%.........
The failure of elements $\infty$ and $\bot$ to satisfy some of these
postulates either invalidates certain standard rules of operation in a CAS
or on the other hand limits the mathematical domain of a CAS. Since the ambition of at least
some CAS is to encompass all of mathematics, such systems must either entirely
march $\bot$ out of mathematics, or ``do something''.
We would like to know that, if we use rules such as $0\times x = 0$,
$x/x=1$ and $x-x=0$, that these statements are true for any possible
legal assignment of values for $x$.
Since these are not true for certain values of $x$, what should we do?
Some of this may be clarified by considering the meaning of equality.
{\bf Or denying the possibility of assigning a value to a name!}
The usual way of providing inference rules for equational logic
includes the following:
Substitution: If $P$ is a theorem then so is $P[x\rightarrow E]$ where
$P[x\rightarrow E]$ means textual substition of expression $E$ for
variable $x$ in $P$.
If $P=Q$ that is, $P$ and $Q$ are of the same type and are equal, then
$E[x\rightarrow P] = E[x\rightarrow Q] $
If $P=Q$ and $Q=R$ then $P=R$.
Some CAS essentially say that two expressions are ``='' under much
looser conditions. For example
\begin{enumerate}
\item[a.] They are equal if, for {\em nearly every} substitution of
values into the variables, the expressions evaluate to equal
expressions. Yes, this ``nearly every'' sounds just wrong, but the
argument is made by the authors of Mathematica that a ``generic
solution'' is acceptable if all exceptions to that solution are in a
space of lower dimensionality. Example: $x/x=1$ even if it is perhaps
false if $x=0$ (That's only one place)
or if $x$ is an interval containing zero (That includes many intervals). Another
example: For all pairs of real numbers $\langle x,y\rangle$ in ${ R}^2$,
$(x^2-y^2)/(x-y)= x+y$ except along a particular line, namely
$x-y=0$). Such an argument defeats the premise of high confidence in
the capability of a CAS program to produce ``guaranteed'' answers.
In particular, it is a hazard if one proposes to use a CAS
prove other programs correct\footnote{We understand why some people
advise against learning
mathematics from a physicist.}.
\item[b.] Two floating-point single- or double-precision
results are equal if they are ``close enough'' (in some relative or
absolute error sense).
\item[c.] Two approximate or imprecise numbers are ``close enough'':
In Mathematica we can construct a very imprecise number $n$,
{\tt n=N[1,1]}, a number 1 with one digit of accuracy. Then
this equation is ``true'': $3+n=4+n$.
\item[d.] Mathematica claims that two values each {\tt Infinity} are equal, yet the
difference of them is {\tt Indeterminate}. This can't be right.
Other examples, by mixing intervals or vectors with scalars, can result
in questionable ``equalities''. [specify?]
\end{enumerate}
How can we reconcile the different behaviors of a CAS resulting from
such designs with the intention that a CAS can ``prove'' mathematical
theorems? If a proof that $f(x)-g(x)$ is zero (or is non-zero) ``for
all $x$'' (where ``prove'' is implemented as simplification), is
contradicted for some $f(c)-g(c)$ at a specific value $c$? Certainly
if $c$ is not in the domain of $x$ we can make some case that this
does not matter. What, though, if it is in the domain?
\subsection{What is equality?}
If we approach this question operationally from the perspective of
a CAS design, the most plausible definition is to say that
two objects are equal if a simplifier program provides the same
canonical result, that is, lexicographically identical, and programmatically
indistinguishable\footnote
{We may (and will) forbid examining the representation at some abstraction
level, otherwise we would be faced with conundrums like declaring as unequal
two integers because they are stored in different memory locations, even if
they had the same value.}, for each side of the relation.
A somewhat weaker demand (for algebras that allow for subtraction) is
that
two objects are equal if their difference simplifies to zero\footnote
{Zero-equivalence algorithms are possibly faster than canonical simplifiers.}.
In Lisp this is {\tt equal}. {Lisp has {\tt eq} to
represent equality as ``same storage location'', as well as several
specialized equalities for numbers ({\tt =}), strings ({\tt string=}),
and combinations ({\tt eql}).} It is possible to have a canonical
expression which is, for example, {\tt 3*inf[3]+14*inf[7]}. which
would be different from {\tt 14*inf[7]+3*inf[3]}; a canonical
simplifier would rearrange one or the other so they would look the
same. Mathematics texts are generally insufficiently pedantic to make
the subtle distinctions necessary for programmers.
We can require that testing for numerical equality, as between
floating-point numbers, use a different predicate. (Advice given to
novice programmers dealing with floating-point computation is to never
test for equality; this advice is overkill: there are
occasions when testing equality does make sense.)
Equality for numbers expressed as intervals has several different
variations---{\em possibly, certainly, impossible}.
\subsection{An alternative approach}
Assume that manipulation in the general simplifier is bound by rules
of computation in a commutative field with identities 0,1 with algebraic
extensions for each independent indeterminate with domain Real or Complex
numbers, and that in this field, $x/x=1$
Any computation that we do must first make sense in this domain. After
the computation is complete, a valuation can be done in some other
domain in which an algebraic expression is associated with an ordered
list of arithmetic commands that can be performed on elements in any
domain, including elements that are infinite, undefined, intervals, or
other objects.
Except for the discomfort caused by $ x^0=1$ and $0^x=0$, this might work.
%***
NEW SECTION 12/2015
Assume we have indexed versions of (positive) infinity oo_i, negative
infinity represented by -oo_i, no assignment to variables
allowed. Complex
infinity - eh, not for now. could be oo*(unit-normal in some
direction) or if no direction specified, just, oh, ooo ?
Anyway, here are some ideas.
oo_1 -oo_2 --> interval[-oo,oo] rather than undefined.
oo_1/oo_2 --> interval [0,oo]
oo_1*oo_2 --> oo_3
oo_1^n --> oo_4 n>1
oo_1^n --> oo_1 n=1
oo_1^n --> 1 n=0
oo_1^n --> 0 n=-1
oo_1^n --> [0, oo] -1 0 n<-1
Mathematica (apparently in the part now known as ``The Wolfram Language'')
has, for several versions, implemented a notion of DirectedInfinity[z]
where the complex number z denotes the direction.
Is |z| normalized?
DI[-1] is printed as -Infinity,
DI[1] as Infinity
DI[] as ComplexInfinity. Should this be DI[0]?
not if z*DI[1] is transformed to DI[z].
Perhaps 0*DI[1] is transformed to indeterminate not DI[0].
Some rules of combination are observed.
DI[1]+DI[-1] is indeterminate, though I think it should be
RealInterval[{-oo,oo}]
consider 8, inverse. $\infty+x= 0$ has no solution. However, if we have a
indexed set of infinities, $\infty_1,~\infty_2,$ etc. then each $\infty_k$ has
an inverse. Real intervals spanning zero arguably do not have inverses
unless the language of intervals is expanded to exterior intervals.
consider 9. If c=inf, then ca=cb becomes inf=inf, or und=und, which does not imply a=b.
If c=und, then ca=cb becomes und=und.
consider 10: und has an unknown sign.
consider 12, 13. the case c=inf or und does not work.
consider 16. und does not have an inverse. inf does not, since inf*0=und, not 1.
Arguably, indexed inf, maybe.
How do intervals fail to satisfy these postulates?
First, what intervals are allowed?
[real,real], [real,inf], [und], [-inf, real], [-inf,inf]
\subsection{Consider how to clarify =}
Two intervals are equal if a simplifier program provides the same
canonical result for each side of the relation. This can happen if
two intervals are each degenerate (one point), and that point is the
same number, or if the two intervals match in upper/lower bounds {\em
and index expression}. Note that in particular it is insufficient
to have equal endpoints alone, although it {\em is} sufficient to
have equal index expression!
consider 8, inverse. $[-1,1] +x= 0$ has no solution. However, if we have
a indexed interval, $[-1,1,p]$, then $[1,-1,-p]$ may work as
an inverse.
consider 9. If c=inf, then ca=cb becomes inf=inf, or und=und, which
does not imply a=b. If c=und, then ca=cb becomes und=und.
%........end insert on algebra
A list of a few plausible rules shows that intervals and the
pseudo-values of $\infty$ and $\bot$ relate to each other. We are here
concerned only with {\em real} values and pseudo-values, and we may
need to have separate rules for {\em complex} $\bot$. This might occur
if one computes $\sqrt{[-4,-1]}$ which has no real value, and could
arguably be $\bot$ but perhaps should be $[1,2]i$ in the context of a
CAS which knows about complex numbers.
\noindent
What is $\sin(\bot)$? [-1,1] if $\bot$ is {\em real}, otherwise $\bot$.\\
What is $\exp(\bot)$? Perhaps $[0,\infty]$ if $\bot$ is {\em real}?\\
Is $[-\infty,\infty]$ the same as real $\bot$?\\
Is $[\bot,\bot]$ meaningful? Is it $\bot$?\\
In general, a plausible result for any function or operation upon $\bot$ is $\bot$.
\medskip
A CAS designer has a choice to make: are these activities to be confined to the
numeric type hierarchy (an uncomfortable choice, I think, in which
case there may be both single- and double-float infinities, and the
choice of projective or affine models of the real line are relegated
to the numeric system), or are such pseudo-values used as an overlay
of a separate algebraic system that coexists more with the language of
variables and constants (like $\pi$ and $e$). The latter seems more
plausible in context of a CAS.
An extended discussion of the consequences of choices to make seems to require
some subjectivity and consideration of applications, and is beyond the scope of the
present paper; a doctrinaire approach is likely to annoy at least some
CAS users. This paper has already exceeded our target length and so
we defer further discussion to detailed program notes
on a particular implementation, described briefly below.
% [[ much more could be added here]]
\subsection{Big-Oh notation}
We hesitate to step into this muck, but one might be tempted to relate
this material to Big-Oh notation in asymptotic analysis. In particular
$O(x)-O(x)$ is $O(x)$, but $O(x) \times O(x)$ is
$O(x^2)$. The common notation is further hindered by the misuse of
the symbol ``='' which is usually considered an equivalence relation.
Yet $T(n)=O(n^2)$ does not mean $O(n^2)=T(n)$; an expression that is
not acceptable. Proposals to remedy this overloading of ``=''
by using curly inequalities or set notation have not driven out the
traditional notation. It would be curious if implementation in a CAS
can provide sufficient motivation to generate an adherence to
notational rigor on mathematicians.
%We now step out of this muck.
%\subsection {Notes on Mathematica}
%[[ particular +/- features from WRI]].
Mathematica has Indeterminate, DirectedInfinity[], and DirectedInfinity[z] for an infinite quantity that is a positive real multiple of the complex number z. Among the rules, 0*Indeterminate is Indeterminate.
%In the typical CAS, including Mathematica, 0*f[x] is simplified to zero; this is a mistake, in the sense that it is false if we later learn that f[x]:=Indeterminate.
%Where do we go with this? A calculus of infinities for Limit computations perhaps. A more cautious simplifier that doesn't say $x/x
%\rightarrow 1$ but $x/x \rightarrow 1 ~~(x\ne 0)$?
%It would be nice to know the rules being used, instead of having to guess them; the treatment has evolved over various versions of Mathematica, and may continue to do so.
%(We could easily omit this section.)
{\section{Implementation notes}
For an arithmetic package in which there are typed objects that look
explicitly like tagged numeric intervals, it is not difficult in
principle to attach methods for all the component-pair combinations
that may be of interest. For example, a two-argument ``+'' in which
each of the arguments is a scalar machine double-float can be
translated fairly easily into an operation. Running through the
plethora of additional operations, e.g. ``+'' of a interval with an
exact rational number, or a symbolic expression like ``f(x+1)'' can be
tiresome, and even with the best of intentions and substantial
effort, inheritance hierarchies can make only some of the right
decision automatically. Indeed, in some cases humans can easily disagree as to
the ``correct'' coercion. Additionally, some unambiguous
operations on ordinary scalars may become unclear with intervals.
For example, there are many more
ways of ``comparing'' intervals: which may involve inclusion, overlap, or
disjointedness. These can be phrased as ``certainly'' or ``possibly''
as in ``certainly greater than''. etc.
Using the Common Lisp Object System (CLOS), we have built a prototype
generic arithmetic system. We define methods for intervals using programs
with headers that look like this:
\begin{verbatim}
(defmethod two-arg-+ ((u Ninterval)(v Ninterval)) ....)
\end{verbatim}
And so forth for all distinct argument type pairs. The traditional
``+'' of any number of arguments is decomposed into adding two arguments
at a time.
That is we can define how to combine any interval
with a complex number, a matrix, or a symbolic variable (or alternatively
refuse to do so). In an object-oriented hierarchy, it
is possible to place intervals in a lattice with respect to some
operations and use inheritance for method resolution, but in some
sense a response like ``{\tt no method available for two-arg-+ with
types Interval and Matrix}'' will wait in the background for
unimplemented (or impossible) combinations. The structure for extending the ordinary
Lisp arithmetic to intervals (namely by overloading existing n-ary
operations like ``+'' and ``*'' by using two-argument versions of
these operations) is straightforward and can be done partly
at compile-time through
macro-expansion\footnote {see generic arithmetic \cite{generic}}
As indicated earlier, using CLOS for keeping track of types does not resolve how to
integrate interval ``stuff'' into a CAS, exactly. This is because we
must patch an existing system, say Maxima, so that it will not apply
identities prematurely that are incorrect for intervals. To
illustrate how pervasive such errors can be, note that {RealInterval} were
introduced in Mathematica in version 3, and even after many years, in
version 5.1, if you type {\tt Interval[{a, b}] - Interval[{a, b}]} you get 0.
which is presumably false. It has been fixed in later versions.
Ideally we have localized a simplification that says $x + (-1)\times x$ =
$0\times x$ = 0.
This would operate only in the case that $x$ is nullable. In
particular $x$ being some version of infinity, or some
as-yet-unspecified interval, makes this simplification false. We must
assert some criterion for $x$ before establishing this exception to
the general simplification. We must also process $x$ to some extent to
see if this criterion holds. It is rather commonplace to try to
``short circuit'' processing by asserting that if any factor in a
product simplifies to zero, the rest of the product need not be
evaluated; this leads to inconsistencies.
A choice remains as to whether we should alter the library of programs
that constitute the ``rational function'' library or not. We could
alternatively argue that these programs are doing arithmetic in a
well-understood field, and we should not use that library if we expect
to be violating the axioms for the ground field of the library, say by
evaluating at an illegitimate point. We could try intercepting any
explicit setting of a rational-variable value to a pseudo-value, but
that would not prevent manipulation of arbitrary symbols, where only
later the value is proffered. Carrying along the assumptions with the
expression (e.g. the answer is ``3'' but only if x is finite), is
a possibility which may be appealing in a theorem proving context, but
contending with the exponential growth in such added conditions on expressions
is daunting. While there does not seem to be a simple solution, at
least recognizing the hazards may make the design of CAS more resilient
to the introduction of such extra values.
While implementation of extra index information on intervals or other
special values solves some of the problems, as indicated above, it also
introduces a complication: Writing index expressions to a file and
reading them back into another running program, or combining them
with expressions from another file requires some way of
distinguishing the index-generated symbols. That is, a data ``dump''
with such indices must re-calibrate them so as not to overlap. This
is similar in nature to using the lisp {\tt gensym} facility in two or more
separate images and avoiding conflicts when combining them).
AHA another idea. Change the CRE (canonical rational expression)
form to allow +-oo, +-0, 0/0. That is, the ratio of polynomials
form allows, with minor extensions, the representation and
somewhat hackish implementation of 1/0, -1/0, 0/1, 0/-1, and for that
matter, other numbers like -6/-1 which is equal to 6 but is more
like limit(x,x,6,minus).
}
%This is done in simptimes, I think.}
\section{Acknowledgments}
A version of this paper was linked to a message on the newsgroup
sci.math.symbolic in June, 2006, which provoked a number of useful
comments. Thanks especially to Christopher Creutzig.
It has been revised as recently as October, 2009, and again, July,
2010.
and yet again, September, 2011, March, 2013, Sept, 2013, June, 2014.
As we update this document in June, 2015
there has been an IEEE approved
draft standard \cite{p1788} which
provides a consensus within the interval-arithmetic
community. This generally places the interval arithmetic
facilities as an extension to numerical floating-point computation,
with complications beyond those expressed above. In this paper
we have ignored the detailed implications in required numerical accuracy to
make intervals fully reliable, the constraints on elementary
functions, interchange formats, and numerous other features discussed
in IEEE 1788-2015.
\begin{thebibliography}{99}
\bibitem {berz-taylor} COSY INFINITY system, http://bt.pa.msu.edu/index\_cosy.htm
%{\em Reliable Computing 4:} 83-97.
\bibitem{emiris} Ioannis Emiris, Richard Fateman,
``Towards an Efficient Implementation of Interval Arithmetic,''
Technical Report UCB/CSD-92-93, July, 1992 (UC Berkeley)
{\tt http://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6247.html}
\bibitem{fate} Richared Fateman, ``To Infinity and Beyond'',
{\verb|http://www.cs.berkeley.edu/~fateman/papers/infinity.pdf|}
\bibitem{hansen}
Eldon R. Hansen, ``Sharpness in Interval Computations'',
{\em Reliable Computing} 3: 17-29 (1997).
\bibitem{p1788} IEEE,
{\tt http://standards.ieee.org/findstds/standard/1788-2015.html}
\bibitem {martin}
Ralph Martin, Huahao Shou, Irina Voiculescu, Adrian Bowyer and Guojin
Wang, ``Comparison of interval methods for plotting algebraic curves,''
Computer Aided Geometric Design, Volume 19, Issue 7, July 2002,
Pages 553-587.
\bibitem{generic} Richard Fateman. file directory \verb|http://www.cs.berkeley.edu/~fateman/generic|
%(http://www.sciencedirect.com/science/article/B6TYN-46C8RF8-1/2/cb147638a60fd464363ae5db7f337ede)
%Author Keywords: Subdivision; Interval analysis; Range analysis;
%Algebraic curves
\bibitem {kahan}
W. Kahan,
``How Futile are Mindless Assessments of Roundoff
in Floating-Point Computation?''
\verb|http://http.cs.berkeley.edu/~wkahan/Mindless.pdf|
\bibitem {berz} Hongkun Liang and Mark A. Stadtherr,
``Computation and Application of Taylor Polynomials
with Interval Remainder Bounds,''
{\tt http://citeseer.ist.psu.edu/cache/papers/cs/8022/ }
\bibitem{markst} http://www.nd.edu/~markst/la2000/slides272f.pdf
\bibitem {mpfi}MPFI group (Revol, Rouillier) INRIA
http://www.cs.utep.edu/interval-comp/interval.02/revo.pdf
\bibitem{nedialko}
Nedialko S. Nedialkov, Vladik Kreinovich, Scott A. Starks,
``Interval Arithmetic, Affine Arithmetic, Taylor Series Methods: Why,
What Next?'' (2003)
{\tt http://citeseer.ist.psu.edu/nedialkov03interval.html}
\bibitem {walster}G. W. Walster, ``Sun Forte Fortran,''
http://developers.sun.com/prodtech/cc/products/archive/whitepapers/tech-interval-final.pdf
\end{thebibliography}
\end{document}
\section{where is the payoff?}
\begin{verbatim}
{do I believe this? consider
$x(x+3)$ vs $x^2+3x$
Let x = [-2,1]
x(x+3) is
[-2,1]*[1,4] = [-8,4]
$x^2+3x$ is
[0,4]+[-6,3] = [-6,7], better, tighter..}
let x= [-2,1]
consider $x*(x+2)$ vs $x^2+2*x$.
[-2,1]*[0,3] = [-6,3] [0,4] + [-4,2] = [-4,6]
A tight bound would be [-1,3]. On [-2,1],
$f(w):=w^2+2*w$ has a min at w=-1, where f(-1)=-1
and a max at 1, where f(1)=3. Our simple routine for {\em NOT}
using Horner's rule
will not work for (x+A)*(x+B). We need to make
more use of symbolic info: the polynomial roots.
\section {How far to extend, e.g. beyond univariate polynomials?}
\begin{verbatim}
Ninterval^NInterval -> no change, or change to Sinterval.
We suspect In the following, k is the next unused interval index.
exp(Ninterval(a,b,#n,f)) -> Ninterval(exp(a){round down},exp(b){round up},#k,1)
log(Ninterval(a,b,#n,f)) -> Ninterval(log(a) etc)
sin(Ninterval(a,b)#n,f)) -> find if rel max min occur in [a,b], check ends
Ninterval(low (round down if not -1), high (round up if not +1) #k, 1)
cos (etc)
must fill out the matrix of programs for all binary ops and unary ops
Ninterval OP other
other OP Ninterval
* - * / expt
=, not=, >, >=, <=, <, included, overlap, disjoint.
perhaps consider reordering of operations if A op B op Ninterval?
\end{verbatim}
-------------------------------------old
Here we are principally interested in integrating into a CAS in some
useful fashion, the manipulation of real number intervals, even given
that they general violate rules of the algebraic system that we are
using for ordinary numbers and symbols. Mathematica (now) does
Intervals in a reasonable fashion, in its latest versions. Maple 7
does not. Maple 7 requires the user to invoke evalr (range arithmetic)
on objects that look like INTERVAL(a..b), or evalrC (complex range).
[check out newer Maple -- I only have Maple 7]. Maple seems to do its
calculations in such a way as to make the assumption that ALL
intervals with the same endpoints represent the same constant, and
thus INTERVAL(-1,1)-INTERVAL(-1,1) becomes zero, even before evalr
does its job.
One consequence of dealing with intervals is that it introduces some
variety in our thinking of the types of computational objects that
might occur in mathematical thinking outside the ring/field extensions
that --while well supported in CAS -- may not be sufficiently
expressive. Resolving what to do about intervals may help in thinking
about computing with other pseudo values Infinity - Infinity as
Indeterminate. \\
Indeterminate - Indeterminate as Indeterminate \\
Interestingly, in MMA 5.1 Interval[{a,b}] - Interval[{a,b}] comes out
0. Oh well. It's been fixed (tested in version 10).
......................
{\section{Other extended numbers}
\section {Symbolic endpoints for intervals}
Intervals with one or two symbolic endpoints, Sinterval(low,hi) are generally
not worth attempting to simplify. They usually get unwieldy, fast.
Perhaps sums are worthwhile, but a product of
two Sintervals is rather clumsy, involving min and max.
An explicit product with a numeric constant may be worth
simplifying: 3*Sinterval(low,hi) --> Sinterval(3*low,3*hi).
An explicit product with a symbol x probably cannot be
simplified, at least if x may be an interval, or not.
Mathematica is rather cautious about its equivalent to
Sinterval, and we could poke around and see what rules it uses.
Intervals with two numeric endpoints (real numbers: single or double
floats, rationals, bignums, bigfloats, QD, MPFR). We exclude
Infinities, NaNs, from the possible endpoints. Given the
context of a symbolic system, we would like to deal with UND
(etc) in a symbolic context, not as a sub-case of an interval.
There are, however, things that can be plausibly most easily
expressed as intervals. For example,
numbers (or alternatively, a representation for ``any real number''
could be an interval from $-\infty$ to $infty$.) Or by an abuse that
seems common, but can be very harmful, this notation
could be ``the set of all possible real numbers.''
Among the other plausible interval concepts: an interval containing
no numbers at all, say [+0,-0] or an exterior interval [1,-1]: the complement of
an ordinary interval [-1,1]. This can be done in an affine or a projective
model (refs); it is possible to represent the union or intersection
of intervals, closed or open endpoints, etc.
We expect that both endpoints should be the same type and precision,
and that given the usual usage of intervals these will be
floating-point. However, rational operations on rational intervals
may be useful, especially since they provide exact (free of rounding)
intervals.}
\end{document}
example from Hongkun Liang and Mark Stadtherr, Univ Notre Dame
AIChe Annual Meeting session 272, Nov. 2000,
f(x) = x ln x for X = [0.3, 0.4]
Using interval arithmetic
F(X) = [0.3, 0.4] * [-1.2040,-0.9163] = [-0.482,-0.275]
Using Taylor model (3rd order, x_0 = 0.35)
F(X) = -0:3674 - 0:04982 *(X - x0) + 1.4286 * (X - x0)^2
-1.3605 * (X - x0)^3 + [-1.25 * 10^-4, 4.86 * 10^-5]
= [-0.370,-0.361]
The exact range is [-0.368,-0.361]
try also
f = sin(2x) + sin(3x) + cos(4x) for X = [0:2; 0:5]
............. do the log example............
3rd order
%keepfloat:true; taylor(x*ln(x),x,0.35,4);
% -0.36744-0.04982*(X-0.35)+1.42857*(X-0.35)^2-1.36055*(X-0.35)^3+1.94364*(X-0.35)^4
$$ -0.36744-0.04982\,\left(x-0.35\right)+1.42857\,\left(x-0.35\right)%
^{2}-1.36055\,\left(x-0.35\right)^{3}+1.94364\,\left(x-0.35\right)^{%
4}+\cdots $$
good answer should be about [1/e, 0.3*ln(0.3)] = appx [-.36788,-.36119185].
note that the 4th derivative of $x \ln x$ is $2/x^3$.
% (* 0.35 (log 0.35)) = -0.36743772
Now try evaluating taylor series expanded at ([0.3,0.4]-0.35) = [-0.05,0.05]
to 3rd order. This is
[-0.05,0.05]
$[-0.05,0.05]^2= [0, 0.0025]$
$[-0.05,0.05]^3= [-1.25\times 10^{-4}, 1.25\times 10^{-4}]$
(+ -0.36744 (* -0.04982 0.05) (* 1.42857 0.0025) (* -1.36055 1.25e-4))
%%= -0.36652964 with error term $1/4! \times 2/\(xi)^3 \times(\xi-0.35)^4$
%% what is the largest value of 2/(xi)^3 for xi in [0.3,0.4]? It is at 0.3, and
%% the value is 74.07407407407406d0
%% that times (1/4!*(0.05)^3
%% which looks like , bounding,..
%% 1/4! * 2/(0.3)^3*0.5^4
%% (* 1/12 (expt 0.3 3)(expt 0.5 4)) = 1.4062e-4
%% (* 1/12 (expt 0.4 3)(expt 0.5 4)) = 3.3333e-4
$1/12\times ([0.3,0.4] -0.35])=$ %??
$0.08333 \times [0.027, 0.064] \times [0, 6.25\times 10^{-4}]= 3.333\times 10^{-4}$
%% (* 0.08333 0.064 6.25e-4) = 3.333\times 10^{-6}$
% (* 0.08333 0.027 6.25e-4) = 1.4061937e-6
[-0.36744, -0.3673066]
%%??? not right. slide gives eror as [-1.25e-4, 4.86e-5]
%%total interval is [-0.370, -0.361]
%%% hummmmmm
{--- this next example is not so convincing --
$$f(x)= \sin(x_0)+\cos(x_0)(x-x_0)-sin(x_0)(x-x_0)^2/2+ \cdots$$ or
using a remainder formula,
where $\xi$ is some value between $x$ and $x_0$, and we (somewhat arbitrarily)
cut off the series after one term:
$$f(x)= \sin(x_0)+\cos(\xi)(\xi-x_0) = \sin(x_0)+ R_s$$
And similarly for $g(x)=\cos(x)$,
$$g(x)= \cos(x_0)-\sin(\xi)(\xi-x_0)= \cos(x_0) + R_c$$
If we then $f(x) \times g(x)$ we get $\sin(x_0)\cos(x_0) + R_s \cos(x_0)+R_c\sin(x_0)$
+ additional remainder terms from the product of the Taylor series which have been
truncated.
Now consider $\sin(x)\cos(x)$ evaluated on $x=[-1,1]$. Since $\sin(x)$
is about [-0.84, 0.84], and $\cos(x)$ is about [0.54, 1], the naive interval
product is [-0.84, 0.84]. Next consider the Taylor model, using
$x_0=0.0$ as the midpoint. Conveniently, $x_0=0,~\sin(x_0)=0,~
\cos(x_0)=1$, so in computing $f\times g$ almost all the terms simplify nicely
leaving us only to compute $ R_s =-\xi \cos(\xi)$, knowing only that
$\xi$ is some value in [-1,1]. A naive computation of this gives us
[-1,1]*[0.54,1] or [-1,1], but if we could do something slightly cleverer, maybe
we could get
[-0.54, 0.54]. This would be much better than [-.84, 84], and not much wider
than the tightest answer, [-0.5, 0.5].
What if we evaluated on $x=[-1/2, 1/2]$ instead?
%-1/4,1/4
then $\sin(x)= [-0.48,48]$
% [-0.248, 0.248]
$\cos(x)= [0.87,1]$.
% [0.9699,1]
The naive interval product is $[-0.48, 48]$.
% [-0.248, 0.248]
Then $ R_s =-\xi \cos(\xi)$, knowing only that
$[-1/2,1/2] \times [0.87,1] = [-1/2,1/2]$
% $[-0.25,0.25] \times [0.96,1] = [-1/4,1/4]$
%DO I BELIEVE THIS??
%CHECK IT OUT again.
%%Carrying another term in the Taylor series:
f(x)=sin(x0)+cos(x0)*(x-x0) - sin(z)*(z-x0)^2/2
g(x)=cos(x0)-sin(x0)*(x-x0) - cos(z)*(z-x0)^2/2
Evaluate f*g on X=[-1,1]
x0=0, sin(x0)=0, cos(x0)=1
f(X) = 0 + x -sin(z)*z^2/2
g(X) = 1 + 0 -cos(z)*z^2/2
f*g= x + - x*(cos(z)+sin(z)) *z^2/2
what can we say about this? looks worse!? z^2/2 = [0,1/2]. cos(z)+sin(z) =? [-.3,1.84]
}
..............
note from G W Walster in Forte doc
consider
A= x/(x+y). This can be rewritten as the SUE
B=1/ (1+y/x)
But the domains of A and B are not the same!
A is defined except when x+y=0.
B is defined except when x=0 or x=-y \ne 0
For intervals,
0 \not\in X+Y,
0 \not\in 1+y/x and 0 \not \in x.
What about
C = 1- (1/(1+x/y)) ?
different domain restrictions!
unless A and B return R* (undef.) this is sharp:
B \cap C % intersection.
Walster asks
When is it legitimate to substitute
one interval expression for another? Is substitution possible when the enclosed
expressions have different domains? The analysis of this and other questions is
facilitated by introducing the closed R* -system in which the containment set of
an interval expression is always defined.
................
some local rules. do we assume inf[a] is positive or unsigned?
This is just writing off the top of my head, and must be checked.
ADDITION
IMPORTANT: assume non-unique inf[a], und[a]
UND RULES
n*und[a]+m*und[a] = (n+m)*und[a] for n,m numbers
n*und[a]+m*und[a] = 0 if n+m=0
in particular, und[a]-und[a] = 0
n*und[a]+ {anything else} = no change
INF RULES
n*inf[a]+m*inf[a]= (m+n)*inf[a]
n*inf[a]+m*inf[a]= 0 if (m+n) = 0
in particular inf[a]-inf[a] =0
n*inf[a]+m*inf[b] = no change?
Assume Q does not involve inf[a].
n*inf[a]+Q = no change?
IMPORTANT: assume UNIQUE inf[a], und[a]
und[a]*+{any} = und[a]
n*inf[a] = inf[a] for n nonzero ; re-use symbol
= und[new] for n zero
inf[a]+inf[b] = und[new]
inf[a]*inf[b] = inf[new] or nochange.
MULTIPLICATION
assume non-unique inf[a]
und[a]*und[a] = no change
und[a]/und[a] = 1 (or und???)
und[a]^n*und[a]^m = und[a]^(n+m) for n,m numbers?
inf[a]^n*inf[a]^m = inf[a]^(n+m) for n,m numbers?
special case und[a]^0? Should it be 1??
special case inf[a]^0? Should it be 1??
und[a]^n = und[new] . what if n=0? plausibly 1??
und[a]* {anything} = und[a]];reuse
inf[a]*inf[b] = inf[c]
inf[a]^n*inf[b]^m = no change
n*inf[a]/m*inf[a]= n/m
n*inf[a]/m*inf[a]= und if m=0 . what if n=0?
n*inf[a]+m*inf[b] = no change?
inf[a]^n*B with n negative, B free of inf[..] or und[..] --> 0
Assume Q does not involve inf.
n*inf[a]+Q = no change?
If inf[a] is unique (no possible collapse with other occurrences of inf[a])
If inf[b] is unique (no possible collapse with other occurrences of inf[b])
n*inf[a] = inf[a] for n nonzero
= und[new] for n zero
inf[a]+inf[b] = und
inf[a]*inf[b] = inf[c] ? or inf[a]? or und?
DIVISION RULES
There are no division rules because $A/B$ is the same as $AB^{-1}$. Look at rules for powers.
SUE for cubics..
consider
3 - 8/(3*Sqrt[3]) + (5 - 2*Sqrt[3])*x^2 + (2/Sqrt[3] + x)^3
which is the same as
3 + 4*x + 5*x^2 + x^3
But has fewer occurrences of x.
In general, for a cubic, eliminate one coefficient ..
$(a+x)^3+b x^2+d=x^3+r x^2+t x+u$
by $\left\{b\to r-\sqrt{3} \sqrt{t},d\to \frac{1}{9} \left(9 u-\sqrt{3}
t^{3/2}\right),a\to \frac{\sqrt{t}}{\sqrt{3}}\right\}$
Examples?
{b -> r - Sqrt[3]*Sqrt[t], d -> (-(Sqrt[3]*t^(3/2)) + 9*u)/9,
a -> Sqrt[t]/Sqrt[3]}
Motzkin polynomials
according to Demmel etal
m4(x,y,z):= z^6+x^2*y^2*(x^2+y^2-4*z^2) has no accurate polynomial
evaluation algorithm.
m3(x,y,z):= z^6+x^2*y^2*(x^2+y^2-3*z^2) does.
That is, there are places for x such that for m4,
computation of p( x,delta) is not zero but p(x)
mathematically is zero.
note x+y+z is not accurately evaluable.
.....
appendix: what could we implement
1. No constant-arithmetic expression is allowed to have, as a subexpression,
oo, -oo, indefinite, empty (the interval [oo,-oo])
It IS possible to have expressions that look like
constant * symbolic expression. e.g. oo*x. or oo+x.
If we represent oo by 1/0 then oo*x prints in Maxima somewhat oddly
as x/0
How to remove the extraneous constant-combinations?
Assuming exactly one occurrence of oo:
(oo)^n for explicit n>0 --> oo
(oo)^n for explicit n<0 --> 0
(oo)^n for not explicit n --> indefinite (arguably could be interval [0,oo]?)
oo+ anything --> oo
-oo+anything --> -oo
sin(oo) --> indefinite
exp(oo) --> oo
exp(-oo) --> 0
atan(oo) --> %pi/2
etc see IEEE 754 literature (or Mathematica, copied from it)
is equal(inf,inf) --> false or unevaluate?
is notequal(inf,inf) --> true
oo*n for explicit non-zero --> oo*sign(n)
oo*n for not explicit n, --> indefinite. could be -oo,oo, or
..mysterious 0*inf
Is it always true that for all f(ind) we should return ind...
sin(ind) --> ind
exp(ind) --> ind
atan(ind) --> ind
sin(inf) --> ind (mathematica says interval [-1,1]
0* ind --> ind
0+ ind --> ind
etc see IEEE 754 literature for NaN
If we encounter subst(a for b in c),
then there are 3 possibilities:
[1] a is oo, -oo, ind
if a "contains" one of these it has to first be reduced to one of
these. Or maybe an interval.
[2] a is an ordinary expression and substituting a for b
produces only an ordinary expression.
[3] when a is substituted, it produces an ind or oo as in
subst(3, b, 1/(b-3)) or subst(3,b,exp(%i*%pi/(b-3))
...
We can try to convert cases [1] and [2] to a limit calculation.
Similarly, an evaluation or a binding could be converted to
a limit calculation as
block([a:oo], x:f(a)) means x:limit(f(z),z,a);
ev(f(a),a=oo) is the same situation.
...............................
What do we do with more than one inf? e.g. consider
limit(2/x^2,x,0) / limit(1/z^2,z,0)
which is oo/oo, probably ind but some might argue is 2.
/*infinity simplification */
load(partition)$
kill(infpat,minfpat,pluspat,timespat,
isinf,isminf,isindef,isplus,istimes,isexpt,
isext, extcount,ext_count,
make_inf,make_indef,
separate_ext_plus,separate_ext_times,
resimp_plus, resimp_times, resimp_expt,
r1,rp,rt,re)$
matchdeclare(infpat, isinf)$
matchdeclare(minfpat,isminf)$
matchdeclare(pluspat,isplus)$
matchdeclare(timespat,istimes)$
matchdeclare([bpat,epat],any)$
isinf(x):= is (not(atom(x)) and ?eq(inpart(x,0),infinity))$
isminf(x):= is (not(atom(x)) and isinf(-x))$
isindef(x):= is (not(atom(x)) and ?eq(inpart(x,0),indefinite))$
isplus(x) := is (not(atom(x)) and (inpart(x,0)="+"))$
istimes(x):= is (not(atom(x)) and (inpart(x,0)="*"))$
isexpt(x):= is (not(atom(x)) and (inpart(x,0)="^"))$
/* maybe....
istimes(x):= is (not(atom(x)) and (inpart(x,0)="*" /* exclude -1*r ??? */
and (not(inpart(x,1)=-1)) and (not (length(x)=2))))$
*/
isext(x):= is (isinf(x) or isminf(x) or isindef(x))$
ext_count:0$
extcount():=ext_count:ext_count+1$
make_inf(z):=infinity(z, extcount())$
make_indef():=indefinite(0,extcount)$
defrule(r1,sin(infpat),make_indef())$ /* for example*/
separate_ext_plus(z):=block([answer],partition_expression
("+",isext,[],cons,"[",answer,z), answer)$
/* returns a list of 2 lists. The first list contains
all extended objects. The second contains everything else. */
separate_ext_times(z):=block([answer],partition_expression
("*",isext,[],cons,"[",answer,z), answer)$
separate_ext_plus (a+indefinite[7]-infinity[5]+infinity[4]+23);
hasindef(L):= some(isindef,L);
hasinf(L) := some(isinf,L);
hasminf(L) := some(isminf,L);
/*apply1(a+b,rp);
apply1(a+infinity[3]+7,rp); */
defrule(rp,pluspat,resimp_plus(separate_ext_plus(pluspat)));
defrule(rt,timespat,resimp_times(separate_ext_times(timespat)));
defrule(re,bpat^epat,resimp_expt(bpat,epat));
resimp_plus(z):=
block([exts:first(z), /* list of extension objects like inf, indef */
regs:second(z)], /* regular numbers, symbols, expressions*/
/* list the cases
exts, regs result
1. [] [a,b,..,z] a+b+...+z
2. [indef,..] any indef
3. [inf(dir,c1)] [a,b,..,z] inf(dir,cl)
4. [inf(dir1,c1),
inf(dir2,c2)] [a,b,..z] if dir1=dir2 then inf(dir1,c3]
else indef
5. more than 2 infs, if dir1=dir2=... then inf else indef
. that's all, I think.
*/
if exts=[] then apply("+",regs)
else
if hasindef(exts) then make_indef() /* indef + any is indef*/
else
if rest(exts)=[] then /* there is only one extension object */
make_inf(inpart(exts,1,1))
/* slight lie?
inf+X is indef if X might later be assigned -inf */
else
/* there is more than one inf object in the sum.
for example oo_1+oo_2 or oo_1-oo2 */
/* if all infs are same direction, addthem? patch later */
if some(isminf,exts) and some(isinf,exts) then make_indef()
else make_inf(inpart(exts,1,1)))$ /* slight lie? */
resimp_times(z):=
block([exts:first(z), /* list of extension objects like inf, indef */
regs:second(z)], /* regular numbers, symbols, expressions*/
/* list the cases
exts regs result
1. [] [a,..,z] a*...*z
2. [indef, ..] [a,...,z] indef
2. [inf(dir,c1)] [a,...,z] inf(signum(dir*a*...*z], c2) ?
or
inf*[signof constants]* [product of symbols and
expressions]
2. [indef, ..] [a,...,z] indef
2. [indef, ..] [a,...,z] indef
3. [inf(dir1,c1),
inf(dir2.c2] [a,...,z] inf(signum(dir1,dir2,a*...z), c3)
4. etc for more than 2 infs.
5. [inf(dir1,c1),
tiny(dir1,c1)] [a,...,z] a*...*z.
This is like 0*inf fully correlated. Cancelling in pairs.
6. [inf(dir1,c1),
tiny(dir2,c2)] [a,...,z] indef. This is like 0*inf, no correlation
7. [inf(dir1,c1),
tiny(dir2,c2),...] [a,...,z] a*...*z indef
8. [inf,..] 0 indef (won't happen because 0
kills it)
9. [tiny] [a,...,z], 0 only one tiny
10. [tiny,tiny..] [a,...,z], 0 only tinys, no indef or inf.
11. we would never see inf(dir,c1)*tiny(dir,c1) because simplifier would
make that into 1. Or inf(dir,c1)^m*inf(dir,c1)^n -->
inf(dir,c1)^(m-n)
which would be inf (m>n) or 0 (m inf(signum(dir1^e),c2)
3. b^e b is indef or inf, e is normal neg
4. b^e b is indef or inf, e is symbolic
5. b^e b is indef or inf, e is 0
6. b^e b is 0, e is symbolic
7. b^e b is 0, e is inf or indef
8. b^e b is 0, e is neg
9. b^e b is 0, e is 0
10.b^(-1) b is inf(dir,c1) --> tiny(dir,c1) ; really?
nan^nan is invalid op in SBCL
floatinf^floatinf is "can't decode nan or inf" in SBCL
0.0*nan = 0.0*floatinf = 0.0
Some of these we could never see because Maxima would
eliminate them, if the simplifier were run on the expression.
Namely
5. inf^0 is changed to 1. und^0 is 1. ind^0 is 1 inf^0.0 is 1.0
6. 0^x is 0 but see 8.
7. 0^und is error, 0^ind is error, but see 8
8. 0 to neg exp error
9. 0^0 is error
BUT WE CAN CHANGE THESE ERRORS. in simpexpt.lisp file.
(-1/0) ^-1 is -0
(1/0) ^ -1 is 0
*/
error("unimplemented: rsimp_expt",base,expon)$
inf_simp(q):=apply1(q,rt,rp)$
/*bug.. inf_simp(-make_inf(0)) loops*/
/* comments ...
infinity(direction, count). count makes each infinity unique.
direction is unit-normal, e.g. 1, -1, or complex direction
as signum(1+%i) is %i/sqrt(2)+1/sqrt(2). So we can have
infinity in a particular complex direction.
infinity(0,count) is 0*oo, which is therefore indefinite (!).
We still need a notation for complex_infinity "all directions".
This comes up for 1/0. (1/(+0) is infinity(1,count). 1/(-0) is
infinity(-1,count). 1/0 is either/both/? or maybe it is
complex_infinity? Here's a thought:
define complex_infinity by infinity(0.0,count).
There is also in IEEE 754 hardware, a signed float zero, which
might require a slight effort to produce, e.g. in lisp (* -1.0 0.0).
We are in trouble with division by 0, not surprisingly. But
here are some details. If 1/0 is replaced by oo, and 1/oo is zero,
what happens to -oo? That is 1/(1/(-oo)) is changed to oo.
Although 1/(-1/oo)= -1/(1/oo) is -oo.
Here's a fix, and it is similar (I hope) to the signed zeros in IEEE
754. (See Kahan paper -- much ado..) In a computer algebra system we
have any number of symbols we can construct, and are not restricted
to 64-bit objects. We could, if we wish, have +0, -0, and unsigned 0.
1/(unsigned 0) would either be the single oo in projective model, or
indefinite since it could be either -oo or +oo in affine. 1/(+0)
would be infinity(1,count), and 1/(-oo) would be infinity(-1,count).
It would be possible to extend this to z/(-oo)
->infinity(signum(-z),count), allowing z to be an arbitrary complex
number or a symbolic expression.
infinity(w1,c1)*infinity(w2,c2) -> infinity(signum(w1*w2),c3)
what about 1/(-oo)? how to keep track? We still have the option
of signed rational numbers which could be 0/(-1), 0/(+1) though this
does not have unsigned 0.
We could introduce, analogous to infinity(..), infinitesimals.
That looks hard to type or pronounce and looks too much like infinity,
so let's use tiny(). That is tiny(1,count) is very much like positive
zero,
tiny(-1,count) is like -0, and tiny(0,count) is unsigned zero.
Reconsider 1/(1/(-oo)), or 1/(1/infinity(-1,c1)).
Running from the bottom up, this becomes 1/(tiny(-1,c2)) and then
infinity(-1,c3)
or -oo. What we wished.
A zero coming from outside this system is treated as tiny(0,ci).
How about 0/0 or 0^0? can we use any of these directions to help??
tiny(a,c1)^n. different cases, e.g. n>0 is tiny(signum(a^n), c2)
tiny(a,c1)^n. n<0 is infinity(signum(a^n), c2)
tiny(a,c1)^0. more properly tiny(a,c1)^tiny(b,c2)... is indef()
given 0/0, we can parse this more carefully looking at the zeros.
Let us say it is tiny(a,c1)/tiny(b,c1) = tiny(a,c1)*tiny(b,c1)^(-1)
= tiny(a,c1)*infinity(signum(b^(-1)),c3).
we need to know a response to ... What is
tiny(a,c1)*infinity(b,c2)?
interval[0, oo] or interval[-oo,0] depending on signum(a*b).
or indefinite() if we want to keep it simple.
Now what is 0^0 just typed in by the user?
tiny(0,c1)^tiny(0,c2). Oh, indef.
Now what is 0/0 just typed in by the user?
tiny(0,c1)/tiny(0,c2). Oh, indef.
note. in Maxima,
ind == indefinite but bounded.
und == undefined as limit(x*sin(x),x,inf)
inf == real positive infinity
minf == real negative infinity (why not -inf ??)
There is an empty set {} internally ($set simp) )
Is the empty set also a real interval? or rather "not an interval"?
We also inherit from the underlying Lisp in those
that have an accessible IEEE 754 floating-point implementation under
them:
single-float
inf,
-inf,
NaNs [many (2^23 of them), some quiet, some signalling]
0.0e0, -0.0e0
double float
inf,
-inf,
NaNs [many (2^52), some quiet, some signalling]
0.0d0, -0.0d0,
maybe extended (long-double) floats, 80 or 128 bit formats
inf,
-inf,
NaNs [many, some quiet, some signalling]
0.0q0, -0.0q0,
CLISP seems to not permit these values.
ARGG.
One suggestion:
convert all float_infs to infinity(1,c) or if minf, infinity(-1,c)
regardless of precision (short,double,extended, long-double, quad).
This will break any effort to do floating ops in the sequel, but they
will probably break, anyway. Any time meval or friends produces a float inf,
nan,
convert it. About -0.0.. uh, leave it alone?
Q. for Maxima users:
what is the relationship between 1/float_inf and 1/infinity(1,c1) ?
many of these questions for IEEE floats and numerical
values, and more are addressed in
Branch Cuts for Complex Elementary Functions
or
Much Ado About Nothing's Sign Bit
by W. KAHAN
https://people.freebsd.org/~das/kahan86branch.pdf
*/
done;
**