| -----Original Message----- | From: boost-bounces@lists.boost.org | [mailto:boost-bounces@lists.boost.org] On Behalf Of John Maddock | Sent: 06 December 2005 18:40 | To: boost@lists.boost.org | 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 +44 1539 561830, Mobile and SMS text +44 7714 330204 mailto: pbristow@hetp.u-net.com http://www.hetp.u-net.com/index.html http://www.hetp.u-net.com/Paul%20A%20Bristow%20info.html