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