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

From: Paul A Bristow (pbristow_at_[hidden])
Date: 2005-12-07 09:45:55

| -----Original Message-----
| From: boost-bounces_at_[hidden] 
| [mailto:boost-bounces_at_[hidden]] On Behalf Of John Maddock
| Sent: 06 December 2005 18:40
| To: boost_at_[hidden]
| Subject: Re: [boost] [math] floating point classification - 
| testinghelpwanted
| 
| > Then for epsilon  2 ^-105 and 2 ^ -106 are below
| > outAsHex(one + numeric_limits<quad_float>::epsilon());
| > // 1.00000000000000000000000000000002 1 == 0x3ff0000000000000
| > 2.4651903288156619e-032 == 0x3960000000000000   2 ^-105
| > // 1.00000000000000000000000000000001 1==0x3ff0000000000000
| > 1.2325951644078309e-032==0x3950000000000000 2 ^ -106
| >
| > I note that the 2 ^ -106 value is the one I would naively expect - a
| > single least significant bit different.
| 
| Strange!
| 
| The docs do make it clear that this type is rather strange:

OK, at least there is no confusion about that!

| Ah, now hold on: when you print out your results, even if there is a 1 in 
| the last binary digit *that doesn't mean there will be a 1 in the last 
| hexadecimal or decimal digit*.  Think about it, unless the number of
binary 
| bits fit into a whole number of bytes, the last binary 1 may 
| be in a 1, 2, 4 or 8 position in the last byte.

But the two values printed are the same in all threee representations,
quad_float, doubles hi and lo, and dougles as hex.  so how are they failing
the comparison.
I've glanced at the quad_float operator> < and == and they look fine.

My primitive routine, growed like topsy, prints the hi and lo doubles in hex
separately, mapping so I don't quite follow how the hex display can not be a
true representation.
(Actually it would be much simpler to union the double to a _int64?) As you
call nextafter repeatedly, on double or quad_float, the bits roll up from
the back end as I would expect.  (Nasty things happen when x.lo hits
numeric_limit<double>::max, but that's surely not he problem here?)

union dhex
{ // Used by void outAsHex(double)
        double d;
        unsigned short usa[4];  // Suits Intel 8087 FP P J Plauger C
Standard p 67.
};

void outAsHex(double d)
{ // displayAsHex
        dhex dd = {d};  // Initialise double dd.d = 0.0;
        // Assume Intel 8087 FP.
        // Standard C Library P J Plaguer p 67.
        char fill = cout.fill('0'); // Save.
        int fmtflags = cout.flags();
        int precision = cout.precision(2 + numeric_limits<double>::digits *
3010/10000); // max_digits10()
        cout << dec << dd.d << "==0x" << hex << right
                << setw(4) << dd.usa[3] // SCCC CCCC CCCC FFFF  sign &
exponent.
                << setw(4) << dd.usa[2] // CCCC
                << setw(4) << dd.usa[1] // CCCC
                << setw(4) << dd.usa[0] // FFFF
                ;
                cout.flags(fmtflags); // Restore.
                cout.fill(fill);
                cout.precision(precision);
}  // void outAsHex(double d)

void outAsHex(quad_float x)
{ // displayAsHex
        char fill = cout.fill('0');
        int fmtflags = cout.flags();
        int precision = cout.precision(2+numeric_limits<quad_float>::digits
* 3010/10000); // Save.
        cout << right << x << ' '; // 33 decimal digits autofloat format.
        outAsHex(x.hi); cout << ' ';	outAsHex(x.lo); cout << endl; //
max_digits10 & hex formats.
        cout.flags(fmtflags); // Restore.
        cout.fill(fill);
        cout.precision(precision);
} //  outAsHex(quad_float x)

I have checked a bit more, but still baffled -- and getting bored ;-)

BOOST_CHECK_SMALL((one + numeric_limits<quad_float>::epsilon()) - one,
epsilon); // Check epsilon with value 'by definition'
// absolute value of (one + numeric_limits<quad_float>::epsilon()) -
one{0.123259516440783094595582588325435e-31} exceeds
0.123259516440783094595582588325435e-31

Which is still bizarre.

outAsHex(one + numeric_limits<quad_float>::epsilon() - one); //
0.246519032881566189191165176650871e-31 2.4651903288156619e-032 ==
0x3960000000000000 0 == 0x0000000000000000
outAsHex(numeric_limits<quad_float>::epsilon()); //
0.246519032881566189191165176650871e-31
2.4651903288156619e-032==0x3960000000000000 0==0x0000000000000000

BOOST_CHECK((one + numeric_limits<quad_float>::epsilon()) != one); // Check
epsilon with value 'by definition' - pass
BOOST_CHECK((one + numeric_limits<quad_float>::epsilon()) >= one); // Check
epsilon with value 'by definition' - pass
BOOST_CHECK(one <= (one + numeric_limits<quad_float>::epsilon())); // Check
epsilon with value 'by definition' - pass

With the tolerance doubled

BOOST_CHECK_SMALL((one + numeric_limits<quad_float>::epsilon()) - one,
epsilon + epsilon); - pass

I conclude so far that this epsilon is a reasonable approximation to an
il--defined value and will _probably_ prove useful.  I'm going to give up
and try it out.

| I actually tried to write a short program to calculate epsilon with this 
| type, but quad_float exhibits some very strange behaviour: 
<snip>
 | Double weird.

Quadruply weird!
 
A proper 128-bit FP would be MUCH nicer.  (This sort of ill-behaviour is
probably why Keith Briggs, 'whose fault all this is ' ;-) now recommends
Exact Real http://keithbriggs.info/xrc.html) despite its lack of speed. Or
NTL RR where you can have enough decimal places to not worry about a few
dodgy ones.  But then espilon is really meaningless?

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