# Boost :

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
| > // 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)

{ // displayAsHex
char fill = cout.fill('0');
int fmtflags = cout.flags();
* 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);

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

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
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

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.

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 +44 1539 561830, Mobile and SMS text +44 7714 330204
mailto: pbristow_at_[hidden]  http://www.hetp.u-net.com/index.html
http://www.hetp.u-net.com/Paul%20A%20Bristow%20info.html
```