But this looks more sensible, but final eps shown is still zero.
Note lo double has become denormalised.
So what is the right value of epsilon - or rather the 'useful' value of eps?
I haven't looked at NTL, but I have been playing with the darwin long double source, and have confirmed that you can can create numbers that have a non-contiguous series of digits. i.e. rather than R = x*2^y we have R = x1*2^y1 + x2*2^y2 where y2 < y1 - std::numeric_limits<double>::digits. If the digits were contiguous then the low part would be normalised so that y2 == y1 - std::numeric_limits<double>::digits in all cases, and the type would behave as a normal floating point number. But as it is, you get the behaviour that Paul's observed: there's no way to calculate epsilon, and maybe any value you calculate is meaningless anyway. I guess if you take epsilon as a measure of error, then the worst case, is also the usual definition of: epsilon = 2^(1-D) If there are D digits in the significand. But the best case, is the same as the formal definition (the difference between one and the next representable number), and gives: epsilon = numeric_limts<double>::denorm_min() So take your pick :-) BTW, I think I can see why they did things this way: if long double is used to temporarily increase the precision of some arithmetic operation, then potentially the format used preserves additional information: particularly for additions and subtractions, it *may* avoid catestrophic cancellation errors. The trouble is you can't easily reason about when it will do it's magic and when it won't. You can't even treat it as a number with twice the precision of a double and ignore the occational extra precision: I suspect (but can't prove)that it leads to the same problems you get with 10-byte long doubles and double-rounding and other issues. It's counter intuitive, but it's well known that there are some operations for which a little bit of extra precision can actually give you the wrong answer (see http://docs.sun.com/source/806-3568/ncg_goldberg.html#3098). Just thought I'd confuse you all some more :-) John.