| Title: | Mathematics at DEC |
| Moderator: | RUSURE::EDP |
| Created: | Mon Feb 03 1986 |
| Last Modified: | Fri Jun 06 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 2083 |
| Total number of notes: | 14613 |
In this topic I show a bug in VAX LISP's multiprecision
integer arithmetic. It would be interesting if those of you
with access to other such packages could try out these
examples.
I came across the example with 2^107 - 2^54 - 1 when a
background batch process tried to verify that the Mersenne
number 2^107 - 1 is indeed prime. The Lucas-Lehmer test
(see note 2.*) showed that it was composite. I immediately
concluded that there was a bug in VAX LISP, that being more
likely than an amazing new discovery about Mersenne primes
[or a bug in my verification program :^)]. Then I isolated
the example. Remember Eric (VIDEO::) Osman's remarks about
verification and my reply in 420.15-16?
The example is extracted from the VAX LISP group's internal
bug tracking conference:
================================================================================
Note 133.3 Bignum multiplication error. 3 of 4
ZFC::DERAMO "Daniel V. D'Eramo" 35 lines 9-DEC-1987 13:03
-< (* bignum bignum) bug in V2.2 >-
--------------------------------------------------------------------------------
There are more bignum multiplication problems. The one
earlier this year was with BIGNUM * FIXNUM; this is a
BIGNUM * BIGNUM problem.
In the following, I have
107 54
a = 2 - 2 - 1
If you square this in VAX LISP V2.2 by evaluating the form
(* a a), the result is
2 160
a - 2
i.e., the result is too small by the given power of two.
Two ways for VAX LISP to compute the correct square are:
(let ((x (expt 2 107))
(y (- (expt 2 54)))
(z -1))
; 2 2 2 2
; (x + y + z) = x + y + z + 2xy + 2xz + 2yz
(+ (* x x) (* y y) (* z z) (* 2 x y) (* 2 x z) (* 2 y z)))
(let* ((a (- (expt 2 107) (expt 2 54) 1))
(r (random a))
(x (- a r)))
; 2 2 2 2
; a = (x + r) = x + 2xr + r
(+ (* x x) (* 2 x r) (* r r)))
The next reply shows how this looks in VAX LISP V2.2.
Dan
================================================================================
Note 133.4 Bignum multiplication error. 4 of 4
ZFC::DERAMO "Daniel V. D'Eramo" 33 lines 9-DEC-1987 13:03
-< Batch job log file showing .-1 >-
--------------------------------------------------------------------------------
$ lisp
Welcome to VAX LISP, version V2.2
Lisp>
(let ((a (- (expt 2 107) (expt 2 54) 1)))
(* a a))
26328072917139289366971320266403017061299610268758208330933993473
Lisp>
(let ((x (expt 2 107))
(y (- (expt 2 54)))
(z -1))
; 2 2 2 2
; (x + y + z) = x + y + z + 2xy + 2xz + 2yz
(+ (* x x) (* y y) (* z z) (* 2 x y) (* 2 x z) (* 2 y z)))
26328072917139290828472957597305935264984442985041227986866536449
Lisp>
(let* ((a (- (expt 2 107) (expt 2 54) 1))
(r (random a))
(x (- a r)))
; 2 2 2 2
; a = (x + r) = x + 2xr + r
(+ (* x x) (* 2 x r) (* r r)))
26328072917139290828472957597305935264984442985041227986866536449
Lisp>
(let ((bad ***) (good1 **) (good2 *))
(list (= bad good1) (= bad good2) (= good1 good2)
(= bad (- good1 (expt 2 160)))))
(NIL NIL T T)
Lisp>
(exit)
$ exit
================================================================================
The earlier BIGNUM * FIXNUM reference is to a bug we think
is fixed in the latest version -- V2.2 -- of VAX LISP. It
would show up as (((415! * 416) * 417) * 418) not being
equal to (415! * (416 * 417 * 418)). I had found this
earlier this year when comparing two different ways of
computing the factorial function.
Dan
P.S. Related notes are:
DIR/TITLE=multiprecision
26 HARE::STAN 2-FEB-1984 2 Brent Multiprecision Package
28 AURORA::HALLYB 4-FEB-1984 0 Integer Multiprecision Package
181 HARE::STAN 18-NOV-1984 0 More Multiprecision
| T.R | Title | User | Personal Name | Date | Lines |
|---|---|---|---|---|---|
| 799.1 | IXP gets the right answer | SQM::HALLYB | Khan or bust! | Thu Dec 10 1987 12:03 | 1 |
| 799.2 | What's IXP? | MATRIX::ROTH | May you live in interesting times | Thu Dec 10 1987 12:50 | 10 |
< Note 799.1 by SQM::HALLYB "Khan or bust!" >
-< IXP gets the right answer >-
What's IXP?
I have MP, but didn't try it. I did try LISP and reproduced the
error...
- Jim
| |||||
| 799.3 | Request for Test Problems | 38464::DERAMO | Now let's see, 1 + 1 = ... | Thu Dec 10 1987 12:53 | 6 |
Does anyone know of a set of test problems with known
answers for debugging bignum arithmetic routines? Should
the test be run before or after the formal program
correctness proof?
Dan
| |||||
| 799.4 | Aha! A pattern emerges. | ZFC::DERAMO | Avid NOTES reader | Thu Dec 10 1987 15:30 | 53 |
; Compile this VAX LISP file and call the MATH function to
; generate 2236 examples of the bug in .0. The search is over
; numbers of the form 2^i +/- 2^j +/- 1. I don't recall seeing
; incorrect squares when both signs were positive. The computed
; squares [when incorrect] were always too small, almost always
; by an integral power of 2^32.
; Without reading the code, can you now determine the bug? :-)
(defun test (x y z i j)
(declare (integer x y) (fixnum z i j))
(let* ((a (+ x y z))
(e (- (the integer (* a a))
(the integer (+ (the integer (* x x))
(the integer (* y y))
(the fixnum (* z z))
(the integer (* 2 x y))
(the integer (* 2 x z))
(the integer (* 2 y z)))))))
(declare (integer a e))
; a = x + y + z = 2^j +/- 2^j +/ 1
; e = a^2 as computed by VAX LISP - the true a^2
(unless (= e 0)
(format t "2^~3D ~A 2^~2D ~A 1 squared is too ~A"
i
(if (> y 0) '+ '-)
j
(if (> z 0) '+ '-)
(if (> e 0) "large" "small"))
(let ((n (truncate (+ 0.5l0 (/ (log (float (abs e) 0.0l0))
(log 2.0l0))))))
(declare (fixnum n))
(if (= (abs e) (expt 2 n))
(format t " by 2^~3D.~%" n)
(format t ".~%"))))))
(defun math ()
(do ((i 75 (1+ i)))
((> i 124))
(declare (fixnum i))
(let ((x (expt 2 i)))
(declare (integer x))
(do ((j 36 (1+ j)))
((> j 70))
(declare (fixnum j))
(let ((y (expt 2 j))
(signs '(-1 1)))
(declare (integer y))
(dolist (sign1 signs)
(declare (fixnum sign1))
(dolist (sign2 signs)
(declare (fixnum sign2))
(test x (the integer (* y sign1)) sign2 i j))))))))
| |||||
| 799.5 | Abridged list of examples from .-1 | ZFC::DERAMO | Avid NOTES reader | Thu Dec 10 1987 15:34 | 39 |
2^ 75 + 2^43 - 1 squared is too small by 2^ 96.
2^ 75 - 2^64 - 1 squared is too small by 2^128.
2^ 75 + 2^64 - 1 squared is too small by 2^128.
2^ 75 - 2^65 - 1 squared is too small by 2^128.
2^ 75 + 2^65 - 1 squared is too small by 2^128.
2^ 75 - 2^66 - 1 squared is too small by 2^128.
2^ 75 + 2^66 - 1 squared is too small by 2^128.
2^ 75 - 2^67 - 1 squared is too small by 2^128.
2^ 75 + 2^67 - 1 squared is too small by 2^128.
2^ 75 - 2^68 - 1 squared is too small by 2^128.
2^ 75 + 2^68 - 1 squared is too small by 2^128.
2^ 75 - 2^69 - 1 squared is too small by 2^128.
2^ 75 + 2^69 - 1 squared is too small by 2^128.
2^ 75 - 2^70 - 1 squared is too small by 2^128.
2^ 75 + 2^70 - 1 squared is too small by 2^128.
2^ 76 + 2^44 - 1 squared is too small by 2^ 96.
2^ 76 - 2^64 - 1 squared is too small by 2^128.
2^ 76 + 2^64 - 1 squared is too small by 2^128.
.
. [2200 lines deleted]
.
2^124 + 2^57 - 1 squared is too small by 2^160.
2^124 - 2^58 - 1 squared is too small by 2^160.
2^124 - 2^58 + 1 squared is too small by 2^ 96.
2^124 + 2^58 - 1 squared is too small by 2^160.
2^124 - 2^59 - 1 squared is too small by 2^160.
2^124 - 2^59 + 1 squared is too small by 2^ 96.
2^124 + 2^59 - 1 squared is too small by 2^160.
2^124 - 2^60 - 1 squared is too small by 2^160.
2^124 - 2^60 + 1 squared is too small by 2^ 96.
2^124 + 2^60 - 1 squared is too small by 2^160.
2^124 - 2^61 - 1 squared is too small by 2^160.
2^124 - 2^61 + 1 squared is too small.
2^124 + 2^61 - 1 squared is too small by 2^160.
2^124 - 2^62 - 1 squared is too small by 2^160.
2^124 - 2^62 + 1 squared is too small.
2^124 + 2^62 - 1 squared is too small by 2^160.
2^124 - 2^63 - 1 squared is too small by 2^160.
2^124 - 2^63 + 1 squared is too small.
| |||||