Dear Prof. Fateman,
my background is computational physics and there large exact powers of two are, well, rare. I understand your arguments, but I don't understand what application would justify the extra effort of supporting such high precision range reduction. Let's look at the case mentioned below of 2^120. Of cause you can represent this number exactly in double precision and you can calculate the sin to 53 bits. But an increment of 1LSB is more than 10^19 pi. I wonder what applications of calculating the sin of a large numbers exist, which jumps angles in increments of more than 10^19 pi. And even more I wonder if such applications are really so important, that they justify the extra effort, from which all applications would suffer?
Don't get me wrong, I fully support a range reduction with high fixed precision, but in order to process any representable power of two, a variable precision range reduction would be required, and this doesn't make sense in my opinion.
Regarding the trap: this gives us a way to know which results are exact and which not, but this doesn't help to resolve the computational effort problem unless there would be a sin_exact function, which does support arbitrary precision range reduction, and another sin function which doesn't. Another problems with traps is, that they are not very feasible with SIMD or parallel processing.
Best regards,
Michael
On 5/11/2012 1:23 AM, Soegtrop, Michael wrote:
Dear Ray,
I slept over your arguments, but I still think they are not right. The problem I see is that floating point numbers can have very large exponents, e.g. more than 1000 binary digits for IEEE double precision. So worst case you need to do a range reduction with a precision of about 1000 bits for a floating point format with 53 bits mantissa.
No, the worst case is that you need more than 1053 bits. How many more depends on how close some multiple of pi is
to an integer that is in range. Say you have x, a number with 1000 bits, an exact power of two, which agrees with
a particular integer multiple of pi to 1018 bits. To do range reduction you need 1053 +18 bits of pi.
The effort for doing the range reduction, essentially a modulus operation, is quadratic (well there are asymptotically faster methods, but I think with SIMD and 1000 bits, the quadratic method is still best), so we have worst case ~400 times the effort, if we assume that numbers are precise.
If you assume the numbers are not precise, why not just discard them and return NaN?
The question is, if this effort is worthwhile. I think no, because I think it is very unlikely, that the result of some computation has a higher precision than the mantissa.
In other words, you are saying all powers of two beyond 53 are not representable, when they are.
The use case of literal numbers is a bit artificial and not really what floating point numbers are made for. One argument for this is that it would be easy to detect if results are exact and add a corresponding flag to floating point numbers.
The IEEE standard has "inexact" traps. Not part of the number but part of the operation.
But this is typically not done, so one can assume that processing of exact numbers is not the intention.
It IS the intention of the trap. It is the unfortunate tendency of language implementors and
hardware implementors to make it difficult or expensive to use this. (By leaving it out of
the language or making it accessible through some excessively complicated mechanism).
So I think, if you need high precision, you should use a high precision floating point implementation.
And if you use the regular routine you should get essentially random?
This is somewhat reminiscent of a method used by some HP calculators
circa 1975? (which allow args up to 10^99) is to implement
a different function from cosine, which looks sort of like cosine but has a period that is not
pi, but an approximation to pi that is something like 12 decimal digits of pi. It has many
of the same properties that cosine has, but has a different period. slightly. This allows
the calculator to maintain some of the trig identities even if the actual results are technically
unrelated to the correct results.
RJF
Best regards,
Michael
From: Raymond Toy [mailto:toy.raymond at gmail.com]
Sent: Thursday, May 10, 2012 6:04 PM
To: Soegtrop, Michael
Cc: Barton Willis; maxima at math.utexas.edu<mailto:maxima at math.utexas.edu>
Subject: Re: cos(1.0d97) with GCL, Clozure, and Julia
On Thu, May 10, 2012 at 8:26 AM, Soegtrop, Michael <michael.soegtrop at intel.com<mailto:michael.soegtrop at intel.com>> wrote:
Dear Barton,
this isn't a problem of the fsin/fcos implementation. For calculating a sine function you need a certain absolute precision. E.g. if you want to have 2 decimals output precision, your number must have an absolute precision of 2 decimals. Floating point numbers have constant relative precision. So the larger the number gets, the worse the absolute precision becomes. If the absolute precision is worse than pi, the result is simply meaningless. Imagine you have the number 1000000000000 with error +-1000. The relative error is quite small (1ppb), but the absolute error is more than 100*pi. The result is meaningless.
I disagree with this analysis. If I want (cos (scale-float 1d0 120)) (from one of Kahan's examples), cos should return the correct value. Sure you can argue that 2^120 (as a float) is some kind of approximation, but, as I've constructed, I want exactly 2^120, which happens to be exactly representable in a double-float. I want cos to return the value as accurately as possible given the input.
Any supposed "error" in a floating-point number is up to the user to decide, The cos function should accept it as an exact number and produce a value based on that.
Ray
_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu<mailto:Maxima at math.utexas.edu>
http://www.math.utexas.edu/mailman/listinfo/maxima
I'm not sure who is taking which position. Here's mine (and I wrote the Maxima cosine bigfloat routine, if that matters...)
Floating point numbers do not have inherent error. Each floating point number is an exact rational number.
2.0^100 is EXACTLY equal to 2^100. It is not 2^100 +- 2^47 with 53 specified digits and the rest "noise". If
you want to represent numbers with error-bands, you can specify the number and the error (2 data items).
Otherwise you are stuck with a base representation (this is done by Mathematica) in which everything is uncertain
and you get a kind of layer upon layer of fuzz in which numbers become anything+-anything. It then becomes
quite tricky to implement a reasonably error system (with 2 data items) because each of the data items
has fuzz!
RJF