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

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