[racket] Error propagation in `flexpt` [was: [racket-bug] all/14862: extflonum returns a different solution to flonum and mixed-integer calculations]

From: Neil Toronto (neil.toronto at gmail.com)
Date: Mon Dec 1 18:17:48 EST 2014

On 12/01/2014 11:46 AM, Mark Lee wrote:
> On 12/01/2014 06:57 AM, Neil Toronto wrote:
>> Moved discussion to the Racket Users mailing list.
>>
>> This isn't an extflonum error. Your extflonum implementation of the
>> fibonacci function is actually more accurate than your flonum
>> implementation.
>>
>> While floating-point operations are fairly accurate, they're not (and
>> can't be) perfect implementations of real functions. Worse, any errors
>> in their inputs may be magnified in their outputs' error, depending on
>> the inputs' actual values.
>>
>> In particular, even if `flexpt` is as accurate as possible (i.e. adding
>> no error beyond rounding to the nearest float), it will magnify error
>> quite a bit when its first or second argument has any error at all and
>> its second argument is large.
>>
>> In this case, you're exponentiating an *approximation* of the golden
>> ratio, which has some error, by a large number. The extra precision in
>> an extflonum makes up for it for the most part, so you get a different
>> answer.
>>
>> If you want to quantify the error in a floating-point function
>> implementation, you can write a separate bigfloat implementation and use
>> as much precision as you like. The default precision of 128 bits is fine
>> in this case.
>>
>>
>> #lang racket
>>
>> (require math/flonum    ; for flround and stuff we'll need later
>>           math/base      ; for phi.0
>>           math/bigfloat)
>>
>> (define (fibonacci-binet-float x)
>>    (flround (fl* (flsqrt #i1/5) (flexpt phi.0 x))))
>>
>> (define (fibonacci-binet-bigfloat x)
>>    (bfround (bf* (bfsqrt (bf 1/5)) (bfexpt phi.bf x))))
>>
>>
>> Here's what I get on your example input:
>>
>>
>>> (fibonacci-binet-float 1474.0)
>> 4.992254605478015e+307
>>
>>> (bigfloat->flonum (fibonacci-binet-bigfloat (bf 1474)))
>> 4.992254605477768e+307
>>
>>
>> Checking it against the exact implementation:
>>
>>
>>> (require math/number-theory)
>>> (fl (fibonacci 1474))
>> 4.992254605477768e+307
>>
>>
>> There are a lot of solutions to this problem. One of my favorites is
>> this one, which uses two floats to represent the golden ratio:
>>
>>
>> (define-values (phi.hi phi.lo) (bigfloat->fl2 phi.bf))
>>
>> (define (fibonacci-binet-float* x)
>>    (flround (fl* (flsqrt #i1/5) (flexpt+ phi.hi phi.lo x))))
>>
>>
>> (See also `make-flexpt`, which is implemented in terms of `flexpt+`.) On
>> your example input, I get
>>
>>
>>> (fibonacci-binet-float* 1474.0)
>> 4.992254605477767e+307
>>
>>
>> This is still one floating-point number away from the closest
>> approximation, but that's about as close as you can get without
>> computing all of it in higher precision.
>>
>> One last thing: floats get incredibly sparse away from zero. In
>> particular, 64-bit floats can exactly represent integers only up to
>> 2^53. (Every float >= 2^53 is even, every float >= 2^54 is divisible by
>> 4, and so on.) Practically, this means it's pretty much impossible for a
>> floating-point Fibonacci implementation to be exact for inputs larger
>> than 74, since (fibonacci 75) > 2^53.
>>
>> For a lot more on how errors propagate through floating-point programs
>> and what to do about it, you can read my article on it, which seems to
>> still be freely available here:
>>
>>
>> http://www.computer.org/cms/Computer.org/ComputingNow/issues/2014/10/mcs2014040080.pdf
>>
>>
>> Neil ⊥
>>
>> On 12/01/2014 12:24 AM, mark at markelee.com wrote:
>>> A new problem report is waiting at
>>>     http://bugs.racket-lang.org/query/?cmd=view&pr=14862
>>>
>>> Reported by Mark Lee for release: 6.1.1
>>>
>>> *** Description:
>>> When running my fibonacci solver, I find that my extfloat code returns
>>> different numbers compared to my flonum code and mixed-integer code.
>>>
>>> *** How to repeat:
>>> #lang racket
>>>
>>> ;; written by Mark Lee
>>>
>>> ;; require some modules
>>> (require racket/flonum)
>>> (require racket/extflonum)
>>>
>>> ;; requires SSE instructions, solve with long double precision
>>> (define (fibonacci-extfloat extfl)
>>>      (if (extfl<= extfl 2.0t0)
>>>          1.0t0
>>>          (extfl+ (fibonacci-extfloat (extfl- extfl 1.0t0))
>>>                  (fibonacci-extfloat (extfl- extfl 2.0t0)))))
>>>
>>> (define (fibonacci-binet-extfloat extfl)
>>>      (extflround (extfl* (extfl/ 1.0t0 (extflsqrt 5.0t0))
>>>                  (extflexpt (extfl/ (extfl+ 1.0t0 (extflsqrt 5.0t0))
>>> 2.0t0) extfl))))
>>>
>>> ;; solve with double precision
>>> (define (fibonacci-float fl)
>>>      (if (fl<= fl 2.0)
>>>          1.0
>>>          (fl+ (fibonacci-float (fl- fl 1.0))
>>>               (fibonacci-float (fl- fl 2.0)))))
>>>
>>> (define (fibonacci-binet-float fl)
>>>      (flround (fl* (fl/ 1.0 (flsqrt 5.0))
>>>                  (flexpt (fl/ (fl+ 1.0 (flsqrt 5.0)) 2.0) fl))))
>>>
>>> ;; test to see if long double and double produce the same result
>>> (= (fl->exact-integer (fibonacci-binet-float (->fl 1474)))
>>>      (extfl->exact-integer (fibonacci-binet-extfloat (->extfl 1474))))
>>>
>>> *** Environment:
>>> Linux x86_64 / Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.21
>>> (KHTML, like Gecko) konqueror/4.14.2 Safari/537.21
>>>
>>
>
> To Neil,
>
> Fascinating stuff, I figured it was an exponentiation issue but wasn't
> sure because my extfloat implementation provided a different answer from
> wolframalpha. Does racket switch precision for operations like * by default?
>
> Regards,
> Mark
>

I'm not sure how to answer the question. When `*` gets a flonum and an 
exact rational, it usually downcasts the exact rational to a flonum and 
then carries out the operation. (The exception is multiplication by 
exact 0, which returns 0. It's weird but makes sense; blame Scheme.)

It could be that your system's 80-bit `pow` isn't correctly rounded. If 
that's the case, there's not really anything we can do to fix it. What 
do you get when you run the following program?


#lang racket

(require racket/extflonum
          math/bigfloat)

(define phi.t 1.6180339887498948482t0)

(real->extfl (bigfloat->real (bfexpt (bf (extfl->exact phi.t))
                                      (bf 1474))))
(extflexpt phi.t 1474.0t0)



The first number should be 1.1163020658834682882t+308. If your system's 
`pow` is correctly rounded, the second number should be the same.

Or it could be that Wolfram Alpha is doing exact real computation.

Neil


Posted on the users mailing list.