From: Ken Dickey
Subject: Numerical gotcha's
Date:
Message-ID: <431@data.UUCP>
I heard an interesting talk recently on interval arithmetic. One of
the messages was about getting bogus answers due to numeric
instabilities but not knowing the numbers were bogus. I am interested
in knowing what results are obtained for the following on various
Scheme & Lisp implementations. Please indicate machine, OS, and what
numerics are available {floats, bignums, rationals}. I will post a
summary [····@data.uucp].
f(x,y) = 333.75 y^6 + x^2 (11 x^2 y^2 - y^6 - 121 y^4 - 2) +
5.5 y^8 + x / (2 y).
for x = 77617 and y = 33096
The Scheme code for the above follows...
;;=====================================================================
;; from talk on interval arithmetic by G. W. Walster, 1990 Nov 19.
(define (f x y)
(+ (* (/ 1335 4) (expt y 6))
(* (expt x 2)
(- (* 11 (expt x 2) (expt y 2))
(expt y 6)
(* 121 (expt y 4))
2)
)
(* (/ 11 2) (expt y 8))
(/ x (* 2 y))
) )
(define (tryit) (f 77617 33096))
;; with flonums is {IBM FORTRAN note wrong sign!}
;; f = 1.17260... single precision
;; f = 1.1726039400531... double precision
;; f = 1.172603940053178... extended precision
;; PC Scheme: f = 1.17260394005318 [Intel 386] {bignums, but not rationals}
;; ELK 1.1: f = 1.18059162071741e+21 [Sun3] ??? {floats only}
;; THE CORRECT ANSWER IS:
;;
;; (tryit) -> -54767/66192 {Gambit 1.4; sun3; rationals+bignums}
;;
;; flonum approx of (- (/ 54767 66192)) is -0.827396059946821 {ELK 1.1}
;;
; ---- e o f ----
;; Thanks,
;; Ken Dickey ····@data.uucp
In article <···@data.UUCP> ····@data.UUCP (Ken Dickey) writes:
I am interested
in knowing what results are obtained for the following on various
Scheme & Lisp implementations.
in t3.1 running on a sun3 under sunos4.1, with a small change because
'-' only takes two arguments, i got
==> (define (f x y)
(+ (* (/ 1335 4) (expt y 6))
(* (expt x 2)
(- (* 11 (expt x 2) (expt y 2))
(+ (expt y 6)
(* 121 (expt y 4))
2)))
(* (/ 11 2) (expt y 8))
(/ x (* 2 y))))
F
==> (TRYIT)
-54767/66192
--
I don't think the stories are "apocryphal". I did it :-) .. ·······@nmsu.edu
In article <···@data.UUCP>, ····@data.UUCP (Ken Dickey) writes:
> ;; ELK 1.1: f = 1.18059162071741e+21 [Sun3] ??? {floats only}
> ;; THE CORRECT ANSWER IS:
> ;; (tryit) -> -54767/66192 {Gambit 1.4; sun3; rationals+bignums}
On an Encore Multimax, T 3.0 gave the correct answer.
However, when the numeric literals were changed to floating point
(e.g. 333.75 instead of (/ 1335 4), 121.0 instead of 121) and the
x, y arguments were written as floats, the answer came out at
roughly -1.1806e21. Since the Sun3 and the Encore both provide
IEEE 754 arithmetic, I imagine other IEEE 754 machines will give
the same answer.
The moral is (if (and (inexact? x) (not (numerical-analyst? *you*)))
(doubt *you* x))
--
I am not now and never have been a member of Mensa. -- Ariadne.
MacScheme 2.0 produces -1.805916207174113e21 which (apparently) is similar to
other IEEE 754 machines. MacScheme has {flonums, bignums} but no rationals.
The manual says, "Flonums are real numbers approximated by a 32-bit floating
point format similar to the IEEE standard." Hmmmm.
Elk 1.2 on a SPARCstation 1 produces the same answer as reported for Elk 1.1
on a Sun3. BTW, Elk has bignums at least.
-- Vladimir
--
==============================================================================
Vladimir G. Ivanovic Sun Microsystems, Inc
(415) 336-2315 2550 Garcia Ave., MTV12-33
········@Sun.COM Mountain View, CA 94043-1100
Disclaimer: I speak only for myself. Your mileage will vary.
==============================================================================
> From: ····@data.UUCP (Ken Dickey)
> Date: 21 Nov 90 17:10:40 GMT
> f(x,y) = 333.75 y^6 + x^2 (11 x^2 y^2 - y^6 - 121 y^4 - 2) +
> 5.5 y^8 + x / (2 y).
> for x = 77617 and y = 33096
(I found your article just this past weekend when searching Google
Groups for interval arithmetic and thought I'd take up your challenge.)
Hmm, the first step for me would be to convert that to LISP
s-expression form, but you seem to have provided the SCHEME translation
already, so there's only a slight change from that for CL:
(defun f (x y)
(+ (* 33375/100 (expt y 6))
(* (expt x 2)
(- (* 11 (expt x 2) (expt y 2))
(expt y 6)
(* 121 (expt y 4))
2))
(* 55/10 (expt y 8))
(/ x (* 2 y))
))
Then just call it (BTW I'm using CMUCL on FreeBSD Unix):
(f 77617 33096) ==> -54767/66192
So what's the difficulty there?
(float * 1d0) ==> -0.8273960599468214d0
;; THE CORRECT ANSWER IS:
;; (tryit) -> -54767/66192 {Gambit 1.4; sun3; rationals+bignums}
;; flonum approx of (- (/ 54767 66192)) is -0.827396059946821 {ELK 1.1}
CMUCL gets the same result, but prints one more significant digit. I
wonder if that extra digit is actually correct. I'll feed that result
into my rational-number formatter and see::
-0.82739605994682136814116509547981629199[0..A]
^^
Yup. ...36...
rounds to ....4
Now to make it more interesting, suppose I'm not allowed to carry exact
rationals through the calculations, but must truncate accuracy to a
limited precision at every point in the calculation. I'll break up the
calculation into individual steps to prove I'm not cheating:
(defun f (x y)
(prog (...)
(setq x2 (expt x 2)) ;All of these are integers in, integers out...
(setq y2 (expt y 2))
(setq y4 (expt y 4))
(setq y6 (expt y 6))
(setq y8 (expt y 8))
(setq c11x2y2 (* 11 x2 y2))
(setq c121y4 (* 121 y4))
(setq bigdif (- c11x2y2 y6 c121y4 2))
(setq c2y (* 2 y))
(setq bigx2 (* x2 bigdif))
(setq c5y8 (* 55/10 y8)) ;Result is integer (bignum)
(setq dxcy (/ x c2y)) ;The only actual fraction which occurs here
(setq c3y6 (* 33375/100 y6)) ;Result is integer (bignum)
(setq mostsum (+ c3y6 bigx2 c5y8)) ;Yields -2, big deal, sigh
(setq bigsum (+ mostsum dxcy)) ;Only other value affected by frac.
(return bigsum)))
(f 77617 33096)
So if we're allowed to have arbitrary bignums (integers), and it's only
rational numbers where we have a limit on precision, there's no change
to the calculation. We get the exactly correct result again,
not worth even running it.
Alternately, to re-interpret the original question: Suppose the
floating point numbers as given are only accurate to the accuracy
shown, rounded to that shown value from some unknown exact value. (But
once we make that concession, we are allowed to do exact
rational-interval arithmetic through the rest of the calculation.)
(setq x2 (expt x 2)) ;All of these are integers in, integers out...
(setq y2 (expt y 2))
(setq y4 (expt y 4))
(setq y6 (expt y 6))
(setq y8 (expt y 8))
(setq c11x2y2 (* 11 x2 y2))
(setq c121y4 (* 121 y4))
(setq bigdif (- c11x2y2 y6 c121y4 2))
(setq c2y (* 2 y))
(setq bigx2 (* x2 bigdif))
;No change to this point. Next the first change:
(setq con5p5
(let ((mid 55/10) (err 5/100)) (list (- mid err) (+ mid err))))
(109/20 111/20)
(setq intvc5y8 (intervals-times con5p5 y8))
(39225688006041672198182894657555398656/5
39945425400647941412828452357694029824/5)
(show-interval *)
7.845137601208334d+36
7.989085080129587d+36
78[4..J]0000000000000000000000000000000000
(setq dxcy (/ x c2y))
77617/66192 ;Same as before of course, but thought I'd show it anyway
(setq con333p75
(let ((mid 33375/100) (err 5/1000)) (list (- mid err) (+ mid err))))
(66749/200 66751/200)
(setq intvc3y6 (intervals-times con333p75 y6))
(10964979499343032646334288396288/25
10965308042976625450200903155712/25)
(show-interval *)
333.745d0
333.755d0
333.74[5..F]
(setq intvmostsum
(intervals-plus
(intervals-plus intvc3y6 bigx2)
intvc5y8))
(-1799343486679944853410296183653957682/25
1799343486679944853410296183653957582/25)
^^^
(setq intvbigsum (intervals-plus intvmostsum dxcy))
(-119102144070318909736934324988422764946519/1654800
119102144070318909736934324988422762208169/1654800)
^^^^^^^
(show-interval *)
-7.19737394671978d+34
7.19737394671978d+34
So given the bounds on the original floating-point constants, we get
extremely wide bounds on possible correct answers. Other packages give
random numbers in that range, whereas interval arithmetic properly
shows the extremely wide range of meaningless garbage out.
Third interpretation: Suppose we assume the constants are accurate all
the way out to the limit of the internal floating point representation,
instead of only as accurate as shown when printed?
(setq f5p5 5.5)
(integer-decode-float f5p5)
11534336
-21
1
(setq imid 11534336)
(setq preint (list (- imid 1/2) (+ imid 1/2)))
(23068671/2 23068673/2)
(setq con5p5 (intervals-divide preint (expt 2 21)))
(23068671/4194304 23068673/4194304)
(show-interval *)
5.499999761581421d0
5.500000238418579d0
5.499999[7..D]
(setq intvc5y8 (intervals-times con5p5 y8))
(7917110997471427464526535127646797564
7917111683866495257675734275403088132)
(show-interval *)
7.917110997471428d+36
7.917111683866495d+36
7917110[9..H]00000000000000000000000000000
(setq dxcy (/ x c2y))
77617/66192 ;Same as before of course, but thought I'd show it anyway
(setq f333p75 333.75)
(integer-decode-float f333p75)
10936320
-15
1
(setq imid 10936320)
(setq preint (list (- imid 1/2) (+ imid 1/2)))
(21872639/2 21872641/2)
(setq con333p75 (intervals-divide preint (expt 2 15)))
(21872639/65536 21872641/65536)
(show-interval *)
333.74998474121094d0
333.75001525878906d0
333.7499[8..C]
(setq intvc3y6 (intervals-times con333p75 y6))
(438605730793681150651170956604
438605770899105173210236705476)
(show-interval *)
4.3860573079368115d+29
4.3860577089910524d+29
4386057[3..8]0000000000000000000000
(setq intvmostsum
(intervals-plus
(intervals-plus intvc3y6 bigx2)
intvc5y8))
(-343197553949286610853411019722
343197553949286610853411019718)
(setq intvbigsum (intervals-plus intvmostsum dxcy))
(-22716932491011179345608982217361007/66192
22716932491011179345608982217251473/66192)
(show-interval *)
-3.4319755394928664d+29
3.4319755394928664d+29
So even if we had full floating-point precision, still the result is
total garbage. Shall I repeat that exercise with double-precision
floats?
(setq f5p5 5.5d0)
(integer-decode-float f5p5)
6192449487634432
-50
1
(setq imid 6192449487634432)
(setq preint (list (- imid 1/2) (+ imid 1/2)))
(12384898975268863/2 12384898975268865/2)
(setq con5p5 (intervals-divide preint (expt 2 50)))
(12384898975268863/2251799813685248
12384898975268865/2251799813685248)
(show-interval *)
5.5d0
5.5d0
5.499999999999999[5..F]
(setq intvc5y8 (intervals-times con5p5 y8))
(1062616696467621908207398403716986070918873023/134217728
1062616696467622079806165352004285857857945665/134217728)
(show-interval *)
7.917111340668959d+36
7.917111340668962d+36
7917111340668960[7..L]00000000000000000000
(setq dxcy (/ x c2y))
77617/66192 ;Same as before of course, but thought I'd show it anyway
(setq f333p75 333.75d0)
(integer-decode-float f333p75)
5871392092323840
-44
1
(setq imid 5871392092323840)
(setq preint (list (- imid 1/2) (+ imid 1/2)))
(11742784184647679/2 11742784184647681/2)
(setq con333p75 (intervals-divide preint (expt 2 44)))
(11742784184647679/35184372088832
11742784184647681/35184372088832)
(show-interval *)
333.75d0
333.75d0
333.7499999999999[7..D]
(setq intvc3y6 (intervals-times con333p75 y6))
(58868667366336962175897158823201458511/134217728
58868667366336972202253164462967895729/134217728)
(show-interval *)
4.386057508463931d+29
4.3860575084639325d+29
4386057508463931[2..A]0000000000000
(setq intvmostsum
(intervals-plus
(intervals-plus intvc3y6 bigx2)
intvc5y8))
(-42899694243660826356810595193/67108864
42899694243660826356542159737/67108864)
(setq intvbigsum (intervals-plus intvmostsum dxcy))
(-177476035086024838637799883019873/277629370368
177476035086024838637340464125537/277629370368)
(show-interval *)
-6.392552590915684d+20
6.392552590915684d+20
Yeah, still gross garbage out given the input tolerance.
The original statement of the problem was ambiguous. Given the nice
form of the numbers, something + 3/4, and something + 1/2,
the two most likely interpretations are:
- They are exactly derived from mathematics (first exact rational calc. above)
- They are crude measurements only accurate to the nearest half or
quarter (even grosser garbage range than the worst I show above)