
Boost : 
From: Paul A Bristow (pbristow_at_[hidden])
Date: 20051204 09:43:49
 Original Message
 From: boostbounces_at_[hidden]
 [mailto:boostbounces_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.4651903288156618919116517665087069677288E32;}; // 2^(1106) == 2^105
0.246519032881566189191165176650871e31 2.4651903288156619e032 ==
0x3960000000000000 0 == 0x0000000000000000
1 + epsilon is
1.00000000000000000000000000000002 1 == 0x3ff0000000000000
2.4651903288156619e032 == 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.2325951644078309e032 == 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 +44 1539 561830, Mobile and SMS text +44 7714 330204 mailto: pbristow_at_[hidden] http://www.hetp.unet.com/index.html http://www.hetp.unet.com/Paul%20A%20Bristow%20info.html
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk