kend@data.UUCP (Ken Dickey) (11/22/90)
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 [kend@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 kend@data.uucpted@nmsu.edu (Ted Dunning) (11/23/90)
In article <431@data.UUCP> kend@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 :-) .. jthomas@nmsu.edubarmar@think.com (Barry Margolin) (11/23/90)
Here are various results on a Symbolics 3600 running Genera 8.0. The "d"
is the double-precision exponent indicator, and "e" is the single-precision
exponent indicator.
Command: (f 77617d0 33096d0)
-1.1805916207174113d21
Command: (f 77617e0 33096e0)
6.338253e29
Command: (f 77617 33096)
-54767/66192
Command: (f 77617 33096e0)
6.338253e29
Command: (f 77617e0 33096)
-6.338253e29
Command: (f 77617d0 33096)
-2.3611832414348226d21
Command: (f 77617d0 33096e0)
1.022335026998684d30
Command: (f 77617e0 33096d0)
-2.0895373854550075d29
--
Barry Margolin, Thinking Machines Corp.
barmar@think.com
{uunet,harvard}!think!barmarok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) (11/23/90)
In article <431@data.UUCP>, kend@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.
vladimir@prosper.EBB.Eng.Sun.COM (Vladimir G. Ivanovic) (11/23/90)
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
vladimir@Sun.COM Mountain View, CA 94043-1100
Disclaimer: I speak only for myself. Your mileage will vary.
==============================================================================kend@data.UUCP (Ken Dickey) (11/27/90)
Thanks to all who responded to my call for numerical results. They came in pretty much as expected. Implementations without rationals (even with bignums) lost. This includes all those based on C. All implementations with rationals apparently also have bignums (not much of a suprise). CL (Common Lisp) implementations are required to have both rationals and bignums. Barry Margolin's Symbolics numbers were especially interesting! It seems that not many numerical analysts are into interval arithmetic. I have not yet heard of a reals implementation based on rationals, despite the advantage of getting error bounds. I have heard that Ray Moore's group at Ohio is getting super-linear performance with some convergence algorithms using interval arithmetic [3+ times speedup upon doubling the number of processors]. Gosh, computers are fun! [I do all my floats on Halloween to frighten the adults 8^]. Thanks again, -Ken kend@data.uucp Implementation result system ============== ====== ====== Scheme->C -1.180591620717411e+21 Sun4/110 SunOS 4.03 C-Scheme 7.0 1.18059162071741e+21 XScheme 1.18059162071741e+21 ELK 1.1 1.18059162071741e+21 Sun3 ELK 1.2 1.18059162071741e+21 SPARCstation 1 MacScheme 2.0 -1.805916207174113e21 T 3.0 approx -1.1806e21 (w reals instead of rationals) Encore Multimax Chez Scheme -54767/66192 SPARC Chez Scheme 3.9g -54767/66192 SPARCstation 1+ SunOS 4.0.3 AllegroCL 3.1.12 -54767/66192 DEC 3100/Ultrix Mac Allegro CL -54767/66192 Mac IIci T -54767/66192 SPARCstation 1 T 3.0 -54767/66192 Encore Multimax T 3.1 -54767/66192 Sun3 SunOS 4.1 Lucid CL -54767/66192 SPARCstation 1 Lucid CL -54767/66192 HP 9000/835 Utah CL -54767/66192 HP 300 BSD 4.3 ================== Here are various results on a Symbolics 3600 running Genera 8.0. The "d" is the double-precision exponent indicator, and "e" is the single-precision exponent indicator. (f 77617d0 33096d0) -1.1805916207174113d21 (f 77617e0 33096e0) 6.338253e29 (f 77617 33096) -54767/66192 (f 77617 33096e0) 6.338253e29 (f 77617e0 33096) -6.338253e29 (f 77617d0 33096) -2.3611832414348226d21 (f 77617d0 33096e0) 1.022335026998684d30 (f 77617e0 33096d0) -2.0895373854550075d29 =================== Gambit 1.5: (digits from Marc Feeley) (floor (* (tryit) (expt 10 100))) -> -8273960599468213681411650954798162919990331157843848199178148416727096930142615421803239062122310854 ;; --- E O F ---