diff --git a/books/bookvol10.1.pamphlet b/books/bookvol10.1.pamphlet index 4a1cd0c..e63cdd5 100644 --- a/books/bookvol10.1.pamphlet +++ b/books/bookvol10.1.pamphlet @@ -3,6 +3,260 @@ \input{bookheader.tex} \mainmatter \setcounter{chapter}{0} % Chapter 1 +\chapter{Interval Arithmetic} +Lambov \cite{Lambov06} defines a set of useful formulas for +computing intervals using the IEEE-754 floating-point standard. + +The first thing to note is that IEEE floating point defaults to +{\bf round-to-nearest}. However, Lambov sets the rounding mode +to {\bf round to $-\infty$}. Computing lower bounds directly uses +the hardware floating point operations but computing upper bounds +he uses the identity +\[\Delta(x) = -\nabla(-x)\] +so that the upper bound of the pair of bounds is always negated. +That is, +\[ x = \left[\underline{x},\overline{x}\right] = +\left<\underline{x},-\overline{x}\right> \] + +Given that convention +\begin{itemize} +\item the sum of $x$ and $y$ is evaluated by +\[\left<\nabla(\underline{x} + \underline{y}), + -\nabla(-\overline{x} - \overline{y})\right> \] +\item changing the sign of an interval $x$ is achieved by swapping +the two bounds, that is $\left<-\overline{x},\underline{x}\right>$ +\item joining two intervals (that is, finding an interval containing +all numbers in both, or finding the minimum of the lower bounds +and the maximum of the higher bounds) is performed as +\[\left\] +\end{itemize} + +Lambov defines operations which, under the given rounding condition, +give the tightest bounds. + +\section{Addition} +\[ x + y = \left[\underline{x}+\underline{y},\overline{x}+\overline{y}\right] +\subseteq \left<\nabla(\underline{x}+\underline{y}) +-\nabla((-\overline{x})+(-\overline{y}))\right>\] +The negated sign of the higher bound ensures the proper direction +of the rounding. + +\section{Sign Change} +\[ -x = \left[-\overline{x},-\underline{x}\right] = +\left<-\overline{x},\underline{x}\right> \] +This is a single swap of the two values. No rounding is performed. + +\section{Subtraction} +\[ x-y = \left[\underline{x}-\overline{y},\overline{x}-\underline{y}\right] +\subseteq \left<\nabla(\underline{x}+(-\overline{y})), +-\nabla((-\overline{x})+\underline{y})\right>\] +Subtraction is implemented as $x+(-y)$. + +\section{Multiplication} +\[ xy = \left[min(\underline{x}\underline{y}, + \underline{x}\overline{y}, + \overline{x}\underline{y}, + \overline{x}\overline{y}), + max(\underline{x}\underline{y}, + \underline{x}\overline{y}, + \overline{x}\underline{y}, + \overline{x}\overline{y})\right] +\] +The rounding steps are part of the operation so all 8 multiplications +are required. Lambov notes that since +\[\Delta(\nabla(r)+\epsilon) \ge \Delta(r) \] +for $\epsilon$ being the smallest representable positive number, one +can do with 4 multiplications at the expense of some accuracy. + +In Lambov's case he makes the observation that +\[ +xy=\left\{ +\begin{array}{l} +\left[min(\underline{x}\underline{y},\overline{x}\underline{y}), + max(\overline{x}\overline{y},\underline{x}\overline{y})\right], + {\rm\ if\ }0 \le \underline{x} \le \overline{x}\\ +\left[min(\underline{x}\overline{y},\overline{x}\underline{y}), + max(\overline{x}\overline{y},\underline{x}\underline{y})\right], + {\rm\ if\ }\underline{x} < 0 \le \overline{x}\\ +\left[min(\underline{x}\overline{y},\overline{x}\overline{y}), + max(\overline{x}\underline{y},\underline{x}\underline{y})\right], + {\rm\ if\ }\underline{x} \le \overline{x} < 0 +\end{array} +\right. +\] +from which he derives the formula actually used +\[xy \subseteq \left +\] +where +\[ +\begin{array}{rcl} +a & = & \left\{ +\begin{array}{rl} +\underline{y} & {\rm if\ } 0 \le \underline{x}\\ +-(-\overline{y}) & {\rm otherwise} +\end{array}\right.\\ +b & = & \left\{ +\begin{array}{rl} +-\underline{y} & {\rm if\ } (-\overline{x}) \le 0\\ +(-\overline{y}) & {\rm otherwise} +\end{array}\right.\\ +c & = & \left\{ +\begin{array}{rl} +-(-\overline{y}) & {\rm if\ } (-\overline{x}) \le 0\\ +\underline{y} & {\rm otherwise} +\end{array}\right.\\ +d & = & \left\{ +\begin{array}{rl} +(-\overline{y}) & {\rm if\ } 0 \le \underline{x}\\ +-\underline{y} & {\rm otherwise} +\end{array}\right. +\end{array} +\] +which computes the rounded results of the original multiplication +formula but achieves better performance. + +\section{Multiplication by a positive number} +If one of the numbers is known to be positive (e.g. a constant) then +\[ {\rm if\ }x > 0 {\rm\ then\ } xy \equiv +\left[min(\underline{x}\underline{y},\overline{x}\underline{y}), + max(\overline{x}\overline{y},\underline{x}\overline{y}) +\right] +\] +This formula is faster than the general multiplication formula. + +\section{Multiplication of Two Positive Numbers} + +If both multiples are positive simply change the sign of the +higher bound on one of the arguments prior to multiplication. +If one of the numbers is a constant this can be arranged to +skip the sign change. + +\section{Division} +Division is an expensive operation. +\[\frac{x}{y} = +\left[ +min\left(\frac{\underline{x}}{\underline{y}}, + \frac{\underline{y}}{\overline{y}}, + \frac{\overline{x}}{\underline{y}}, + \frac{\overline{x}}{\overline{y}}\right), +max\left(\frac{\underline{x}}{\underline{y}}, + \frac{\underline{x}}{\overline{y}}, + \frac{\overline{x}}{\underline{y}}, + \frac{\overline{x}}{\overline{y}}\right) +\right]\] +which is undefined if $0 \in y$. To speed up the computation Lambov +uses the identity +\[ \frac{x}{y} = x\frac{1}{y} \] + +Lambov does a similar analysis to improve the overall efficiency. +\[\frac{x}{y} = \left\{ +\begin{array}{rl} +\left[min(\frac{\underline{x}}{\underline{y}}, + \frac{\underline{x}}{\overline{y}}), + max(\frac{\overline{x}}{\underline{y}}, + \frac{\overline{x}}{\overline{y}})\right], & +{\rm if\ }0 < \underline{y} \le \overline{y}\\ +{\rm exception} & +{\rm if\ }\underline{y} \le 0 \le \overline{y}\\ +\left[min(\frac{\overline{x}}{\underline{y}}, + \frac{\overline{x}}{\overline{y}}), + max(\frac{\underline{x}}{\underline{y}}, + \frac{\underline{x}}{\overline{y}})\right], & +{\rm if\ }\underline{y} \le \overline{y} \le 0 +\end{array} +\right. +\] + +The formula he uses is +\[\frac{x}{y} \subseteq +\left +\] +where +\[ +\begin{array}{rcl} +a & = & \left\{ +\begin{array}{rl} +\underline{x} & {\rm if\ }(-\overline{y}) \le 0\\ +-(-\overline{x}) & {\rm otherwise} +\end{array}\right.\\ +b & = & \left\{ +\begin{array}{rl} +(-\overline{x}) & {\rm if\ }0 \le \underline{y}\\ +-\underline{x} & {\rm otherwise} +\end{array}\right. +\end{array} +\] + +\section{Reciprocal} +\[\frac{1}{x} = \left[\frac{1}{\overline{x}},\frac{1}{\underline{x}}\right] +\subseteq +\left<\nabla\left(\frac{-1}{(-\overline{x})}\right), + \nabla\left(\frac{-1}{\underline{x}}\right) +\right> +\] +which is undefined if $0 \in x$. Lambov implements this by checking for +zero, followed by division of $-1$ by the argument and swapping the +two components. + +\section{Absolute Value} + +\[\abs{x} = +\left[max(\underline{x},-\overline{x},0), + max(-\underline{x},\overline{x})\right] = +\left +\] + +\section{Square} +\[x^2 = \abs{x}\abs{x}\] +using multiplication by positive numbers, mentioned above. + +\section{Square Root} + +\[\sqrt{x} = \left[\sqrt{\underline{x}},\sqrt{\overline{x}}\right]\] +which is defined if $0 \le \underline{x}$ + +Lambov notes that this formula has a rounding issue. He notes that +since +\[\Delta(r) \le -\nabla(-\epsilon - \nabla(r))\] +he uses the formula +\[\sqrt{x} \subseteq +\left\{ +\begin{array}{l} +\left<\nabla\left(\sqrt{\underline{x}}\right), + -\nabla\left(\sqrt{-(-\overline{x})}\right)\right>, +{\rm\ if\ }\nabla\left(\nabla\left(\sqrt{-(-\overline{x})}\right)\right)^2 += -(-\overline{x})\\ +\left<\nabla\left(\sqrt{\underline{x}}\right), + \nabla\left(\nabla\left(-\epsilon-\sqrt{-(-\overline{x})}\right) + \right)\right>, +{\rm\ otherwise} +\end{array} +\right. +\] +where $\epsilon$ is the smallest representable positive number. + +The first branch of this formula is only satisfied if the result of +$\sqrt{-(-\overline{x})}$ is exactly representable, in which case +\[\nabla\left(\sqrt{-(-\overline{x})}\right) = + \nabla\left(\sqrt{-(-\overline{x})}\right) +\] +otherwise the second branch of the formula adjusts the high bound +to the next representable number. If tight bounds are not required +the second branch is always sufficient. + +If the argument is entirely negative, the implementation will raise +an exception. If it contains a negative part, the implementation will +crop it to only its non-negative part to allow that computations +such as $\sqrt{0}$ ca be carried out in exact real arithmetic. + + \chapter{Integration} An {\sl elementary function} \index{elementary function} @@ -8355,6 +8609,10 @@ PhD thesis, MIT, Computer Science, 1984 M. van Hoeij. ``An algorithm for computing an integral basis in an algebraic function field'' {\sl J. Symbolic Computation} 18(4):353-364, October 1994 +\bibitem [Lambov 06]{Lambov06} Lambov, Branimir\\ +``Interval Arithmetic Using SSE-2\\ +in Lecture Notes in Computer Science, Springer ISBN 978-3-540-85520-0 +(2006) pp102-113 \bibitem[Wa03]{Wa03} Watt, Stephen, ``Aldor'', \verb|www.aldor.org| \bibitem[We71]{We71} diff --git a/books/bookvolbib.pamphlet b/books/bookvolbib.pamphlet index 0be4854..2b24ac7 100644 --- a/books/bookvolbib.pamphlet +++ b/books/bookvolbib.pamphlet @@ -3194,5 +3194,12 @@ Greve, David A.; McClurg, Jedidiah R.\\ ``Development of a Translator from LLVM to ACL2''\\ \verb|arxiv.org/pdf/1406.1566| +\subsection{Interval Arithmetic} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\bibitem [Lambov 06]{Lambov06} Lambov, Branimir\\ +``Interval Arithmetic Using SSE-2\\ +in Lecture Notes in Computer Science, Springer ISBN 978-3-540-85520-0 +(2006) pp102-113 + \end{thebibliography} \end{document} diff --git a/changelog b/changelog index fa1c2dc..263703c 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,6 @@ +20140617 tpd src/axiom-website/patches.html 20140617.01.tpd.patch +20140617 tpd books/bookvol10.1 add interval arithmetic +20140617 tpd books/bookvolbib add interval arithmetic 20140616 tpd src/axiom-website/patches.html 20140616.01.tpd.patch 20140616 tpd src/input/Makefile add polygamma 20140616 tpd src/input/polygamma.input add a polygamma example diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index 263a145..b76c633 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -4474,6 +4474,8 @@ buglist bug 7258: acosh(0.0) invalid argument to acosh buglist wish 1011: sum(1/(k+a), k=1..n) by Gosper's method 20140616.01.tpd.patch src/input/polygamma.input add a polygamma example +20140617.01.tpd.patch +books/bookvol10.1, bookvolbib add interval arithmetic