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, m...@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
