Re: [boost] [math] floating point classification - testinghelpwanted

From: Paul A Bristow (pbristow_at_[hidden])
Date: 2005-12-04 09:43:49

| -----Original Message-----
| From: boost-bounces_at_[hidden] 
| [mailto:boost-bounces_at_[hidden]] On Behalf Of John Maddock
| Sent: 03 December 2005 19:09
| To: boost_at_[hidden]
| Subject: Re: [boost] [math] floating point classification - 
| testinghelpwanted
| 
| > ld = hi + low as above is not among the bit patterns he lists as
| > valid. If you apply the definition of epsilon to the "normal" long
| > doubles only, then you get a value like the one Paul 
| Bristow computed
| > for NTL's quad_float.
| 
| I agree: that explanation is quite explicit when it says it 
| follows Kahan's 
| "double double" and the IEEE spec.  So the low part must be 
| normalised so 
| that it's bit's "follow on" from those in the high part.  As 
| it stands, 
| 1+numeric_limits<>::epsilon() should evaluate to 1, but we 
| really need to 
| check this out.  Does anyone of a numerical inclination, want 
| to run some tests on Darwin?

Am I right in concluding that Darwin is same as Keith Briggs doubledouble
and NTL?

http://www.cs.cornell.edu/Info/People/vavasis/qmg2.0/patch2_files/doubledoub
le2.h

and www.shoup.net/ntl/ gives references to doubledouble and quad_float.

For NTL quad_float, using my 'guess' at epsilon.

using NTL::quad_float;
static quad_float epsilon() throw() { return
2.4651903288156618919116517665087069677288E-32;}; // 2^(1-106) == 2^105

0.246519032881566189191165176650871e-31 2.4651903288156619e-032 ==
0x3960000000000000 0 == 0x0000000000000000

1 + epsilon is
1.00000000000000000000000000000002 1 == 0x3ff0000000000000
2.4651903288156619e-032 == 0x3960000000000000

So ls digits is 2 - twice what I would expect.

What I get from my nextafter(1) function is

1.00000000000000000000000000000001 1 == 0x3ff0000000000000
1.2325951644078309e-032 == 0x3950000000000000

Note now get a 1 as the least significant bit.

So, perhaps naively, I suspect my epsilon is twice as big as it should be.

quad_float one;
one = to_quad_float("1.");
BOOST_CHECK(one + numeric_limits<quad_float>::epsilon() == one);

./unitTestFunc2.cpp(556): error in "log10_test_function": check one +
numeric_limits<quad_float>::epsilon() == one failed

as I would expect.

I'd also like to have nextafter (and the other IEEE recommended) for
quadfloat.

My tentative stab at nextafter for quad_float is

quad_float nextafter(const quad_float& x, double y =
numeric_limits<double>::max())
{ // Change by one unit in last place (least significant bit), default
increases.
        // Might check that x (and y) are finite?
        // What about sign of x and y?
        double direction = (x < y) ?  +numeric_limits<double>::max() :
-numeric_limits<double>::max();
        // If x < y then increment by one ulp, else decrement.
        quad_float temp = x;
        if (temp.lo == numeric_limits<double>::max())
        { // nextafter would overflow.
                temp.lo = 0.;
                temp.hi = std::tr1::nextafter(temp.hi, direction);
        }
        else
        { // Normal case.
                temp.lo = std::tr1::nextafter(temp.lo, direction);
        }
        return temp;
} // quad_float nextafter(const quad_float& x)

Paul

-- 
Paul A Bristow
Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB
Phone and SMS text [phone], Mobile and SMS text [phone]
mailto: pbristow_at_[hidden]  http://www.hetp.u-net.com/index.html
http://www.hetp.u-net.com/Paul%20A%20Bristow%20info.html