| 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
![[logo]](/logo.png){kind=logo}