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 floatingpoint 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 floatingpoint 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 (fibonaccibinetfloat x) >> (flround (fl* (flsqrt #i1/5) (flexpt phi.0 x)))) >> >> (define (fibonaccibinetbigfloat x) >> (bfround (bf* (bfsqrt (bf 1/5)) (bfexpt phi.bf x)))) >> >> >> Here's what I get on your example input: >> >> >>> (fibonaccibinetfloat 1474.0) >> 4.992254605478015e+307 >> >>> (bigfloat>flonum (fibonaccibinetbigfloat (bf 1474))) >> 4.992254605477768e+307 >> >> >> Checking it against the exact implementation: >> >> >>> (require math/numbertheory) >>> (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: >> >> >> (definevalues (phi.hi phi.lo) (bigfloat>fl2 phi.bf)) >> >> (define (fibonaccibinetfloat* x) >> (flround (fl* (flsqrt #i1/5) (flexpt+ phi.hi phi.lo x)))) >> >> >> (See also `makeflexpt`, which is implemented in terms of `flexpt+`.) On >> your example input, I get >> >> >>> (fibonaccibinetfloat* 1474.0) >> 4.992254605477767e+307 >> >> >> This is still one floatingpoint 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, 64bit 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 >> floatingpoint Fibonacci implementation to be exact for inputs larger >> than 74, since (fibonacci 75) > 2^53. >> >> For a lot more on how errors propagate through floatingpoint 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.racketlang.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 mixedinteger 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 (fibonacciextfloat extfl) >>> (if (extfl<= extfl 2.0t0) >>> 1.0t0 >>> (extfl+ (fibonacciextfloat (extfl extfl 1.0t0)) >>> (fibonacciextfloat (extfl extfl 2.0t0))))) >>> >>> (define (fibonaccibinetextfloat 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 (fibonaccifloat fl) >>> (if (fl<= fl 2.0) >>> 1.0 >>> (fl+ (fibonaccifloat (fl fl 1.0)) >>> (fibonaccifloat (fl fl 2.0))))) >>> >>> (define (fibonaccibinetfloat 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>exactinteger (fibonaccibinetfloat (>fl 1474))) >>> (extfl>exactinteger (fibonaccibinetextfloat (>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 80bit `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