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

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
>