Re: comparing doubles for equality

From: Richard Weeks <rwe...@xxxxxxxxxx>
Is it possible to "just use" interval arithmetic. That is, is
there a refined and established method that satisfactorily deals
with the intrinsic problems of interval arithmetic, e.g. an
interval that includes zero?

There's nothing special about an interval that includes zero,
except that multiplication of two intervals that *both* include
zero requires four basic multiplications instead of the usual two,
followed by computing MIN or MAX of the appropriate possible lower
and upper bounds respectively. That's a trivial amount of extra
work which hardly ever is needed. Division by an inteval that
includes zero is very simple: Throw an exception! If you divide
by a number whose sign you don't even know for sure, then you know
for sure your algorithm is worthless at that point, and you should
definitely bail out. (By comparison, with ordinary floating-point
computations, most of the time you get a number very close to zero
on either side, essentially random which side you are on, and
division gives a very large quotient, which swamps all the rest of
your valid numbers, giving you totally weird results, with no way
to trace it back to the division-by-effectively-zero that happened
earlier.)

The key thing is that you *must* be careful to use interval
arithmetic all the way from input. For example, if you input an
exact value of 4.7, and your internal representation is an interval
of regular floats, you do *not* simply use floating-point input,
giving a value almost but not exactly 4.7, and then compute from
that exact value. You *must* input the 4.7 as an interval of two
adjacent floating-point values, and then feed that interval through
your calculations.

And if you input a value which is not exactly known, you *must*
give it as an explicit interval somehow. For example, if you
measure somebody's under-the-tongue temperature of 98.9 degrees
fahrenheit, you need to know how accurate that actually is,
probably the range of possibilities is 98.8 to 99.0, so you enter
that interval instead of the exact value of 98.9.

Once you get all your input values properly defined as intervals,
then you simply do all your arithmetic (except exact integers such
as indexes of computing loops and indexes of arrays etc.) using
interval arithmetic. Then your output must also be interval
arithmetic. For example, of your final result is between 16.43 and
16.48, you might display that as 16.4[3..8] to make it easy to
visually see which are significant digits and which are not
precisely known. I like that notation better than the standard
notations of 16.455 +/- 0.025 or [16.43 .. 16.48], but any of them
might be your personal favorite.

Isn't interval arithmetic still a work in progress?

Well yes, in the sense that AFAIK none of the major high-level
programming languages (Common Lisp, Java, C++, Perl, ...)
provide it as part of the built-in library. For each language you
can get it only via some third-party add-on (or write it yourself).
AFAIK there is no ANSI/ISO standard for it.
Before I post, let me do a Google search for that ...

<http://developers.sun.com/sunstudio/products/previous/fortran/features.html>
Fortran 95 seem to have some built-in operators but it's not clear
that it has a complete interval-arithmetic library.
Also:
This release of the C++ compiler provides an interface to the interval
arithmetic library. ...
That seems to be just one vendor's C++ compiler that has that
enhancement. If your program is already in C++, you might consider
switching to Sun's C++ compiler with interval arithmetic and see
how it deals with input and output (from the blurb I snipped, it
seems to handle internal calculations just fine).
Also:
The C++ interval library is now compatible with interval arithmetic,
as implemented in the Fortran compiler. ...

<http://ra.dkuug.dk/jtc1/sc22/open/n2848.htm>
V8 - (Interval Arithmetic Requirement)
J3 has found difficulties with the proposed content and requested
further guidance.
(It sounds like ISO is working toward a standard but doesn't have
it yet, but this if FORTRAN only, not any language you'd want to
use for high-level software.)

<http://www.oonumerics.org/oon/>
* Gaol, C++ library for interval arithmetic. Includes methods for
interval constraint solvers.
Looks like a commercial third-party add-on for C++.
* C-XSC: A C++ Class Library for Extended Scientific Computing (A
C++ interval methods class library)
* LIA,GIA,ICE libraries for interval methods in C++ from Delisoft
Ltd. Includes interval arithmetic, global optimization, and
solving nonlinear equations.
Ditto, ditto.

Google shows several more articles about Fortran 95 with interval arithmetic.

<http://portal.acm.org/citation.cfm?id=227701&coll=portal&dl=ACM>
43 KAHAN, W. M. 1980. Interval arithmetic options in the proposed
IEEE Floating Point Standard. In Proceedings of the International
Symposium on Interval Mathematics. Academic Press, New York, 99-128.
Note that just about all modern programming language standards
include IEEE-standard floating-point numbers, together with some
hacks (definitely not thread safe!!) to set rounding mode to
lowerBound or upperBound or boundNearestZero or
boundFurthestFromZero etc. just before doing a primitive arithmetic
operation, the first two of which can be combined to craft together
your own interval arithmetic package. But note that these rounding
modes do *not* do the "right thing" if you call anything other than
a primitive arithmetic (add/sub/mul/div) operation. In particular
if you try these with trig functions you get a totally trashed
result because *every* single operation within the truncated-Taylor
series approximation suffers the *same* rounding mode you have set
globally. You need to completely re-write all the trig library
functions to set the rounding mode interally to get the correct
result. In summary, IEEE rounding modes is nowhere close to having
such an interval arithmetic package already supplied, nor even
formally specified by IEEE.

<http://oracle-wtf.blogspot.com/2006/02/stop-press-oracle-granted-license-to.html>
Yes, it's official. Oracle has been granted permission to extend
February by 3 days. Shame no-one told the developers responsible for
INTERVAL arithmetic. ...
(OT joke, based on accurately implementing the faulty ANTI/ISO SQL standard)

<http://stommel.tamu.edu/~baum/linuxlist/linuxlist/node24.html>
Aha, a succinct definition of interval arithmetic and explanation
why it's such a good thing:
interval arithmetic
A type of arithmetic defined on sets of intervals rather than
sets of real numbers. The usefulness of this sort of arithmetic
lies in its implementation on computers, e.g. outwardly rounded
interval arithmetic allows rigorous enclosures for the rnages
of operations and functions. This simply means that results in
scientific computations are intervals rather than numbers, a
qualitatively different situation than before since the result
is guaranteed to lie in the calculated interval rather than
probably be located somewhere near the single result. Interval
arithmetic also has some utility in automatic theorem proving.
Next:
Software packages which implement interval arithmetic include:
(list of 8: several C++ libraries, one C library with C++
interface, several Fortran libraries, one Matlib library)

<http://ralyx.inria.fr/2005/Raweb/arenaire/uid91.html>
The challenges encountered when developing the Boost library made
clear how an interval arithmetic library and the C++ language have to
interact. As a consequence, G. Melquiond wrote, in collaboration with
H. Brvnnimann (Polytechnic U. Brooklyn, NY USA) and S. Pion
(Giomitrica team, Sophia Antipolis), a proposal [50] for integrating
interval arithmetic to the C++ standard library which was submitted to
the C++ Standards committee (ISO/IEC JTC1/SC22/WG21).
->
[50]
H. Brvnnimann, G. Melquiond , S. Pion.
A Proposal to add Interval Arithmetic to the C++ Standard
Library, Technical report, C++ standardization committee, 2005,
n^o 5646 http://www.inria.fr/rrrt/rr-5646.html.
-> <http://www.inria.fr/rrrt/rr-5646.html>
Abstract : Interval arithmetic is a basic tool for certified
mathematical computations, it is presented in many references. We
describe here the formal proposal to include interval arithmetic in
the C++ standard library.

In summary, I was mostly correct, no standard yet, only a bunch of
third-party add-ons for C++ and Fortran (and C and Matlib one
each), plus one major ventor (Sun) which has it in their Fortran 95
implementation.

Java or Common Lisp? Roll your own, sigh.

Actually I've been for quite some months proposing getting together
to craft a cross-language specification for interval arithmetic, to
be implemented in both Common Lisp and Java initially. But I
haven't been able to find anybody interested in working with me on
either the draft standard *or* any implementation work. I can't
even find "guinea pigs" to test what I write myself on my own time
without pay, so I've had no incentive to do any further work past
what I did 2004.Apr:
<http://www.rawbw.com/~rem/IntAri/>
Note especially the comparison between IEEE double-precision
floating point and the interval arithmetic I implemented myself:
<http://www.rawbw.com/~rem/IntAri/compareFltIntArith1.txt>
.

Relevant Pages

  • Re: C89, size_t, and long
    ... libraries which have useful stuff). ... So because others here can't think of a good reason to use long instead of size_t and you can't prove there is never such a reason, we must all assume there could be a good reason to not use the mechanisms provided by the standard to deal with a problem? ... If you are claiming there is a reason to avoid the mechanisms the standard provides to allow portability then it is up to *you* to prove your point, not up to others to disprove it. ... but you can hardly blame the consequences on the Standard headers. ...
    (comp.lang.c)
  • Re: C89, size_t, and long
    ... libraries which have useful stuff). ... So because others here can't think of a good reason to use long instead of size_t and you can't prove there is never such a reason, we must all assume there could be a good reason to not use the mechanisms provided by the standard to deal with a problem? ... want to have standard headers included in their own headers which ...
    (comp.lang.c)
  • Re: "Criticism of the C programming language ??????"
    ... "that's not in the standard, ... rewrites of the compiler. ... PH> They got all the compiler vendors to agree to a common ... or graphics libraries or threading libraries, ...
    (comp.lang.c)
  • Re: OT: Requesting C advice
    ... fit in the standard process memory footprint of a VAX at the time. ... in the libraries--and indeed how many libraries there were going to be. ... had (a PDP-11 and a VAX), the linker couldn't handle the library size. ... a nice Morrow Microdecision with two floppy disk drives. ...
    (Fedora)
  • Re: OT: Requesting C advice
    ... C was NEVER a "virtual machine language" system. ... But the Standard uses the term "abstract machine", ... libraries was not pleasant. ... that was a huge bone of contention in the committee meetings. ...
    (Fedora)