1 ;;;; -*- Mode: lisp -*-
3 ;;;; Copyright (c) 2007 Raymond Toy
5 ;;;; Permission is hereby granted, free of charge, to any person
6 ;;;; obtaining a copy of this software and associated documentation
7 ;;;; files (the "Software"), to deal in the Software without
8 ;;;; restriction, including without limitation the rights to use,
9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10 ;;;; copies of the Software, and to permit persons to whom the
11 ;;;; Software is furnished to do so, subject to the following
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
28 ;;; double-double float routines needed for quad-double.
30 ;;; Not needed for CMUCL.
32 ;;; These routines were taken directly from CMUCL.
34 (declaim (inline quick-two-sum
))
35 (defun quick-two-sum (a b
)
36 "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
37 (declare (double-float a b
))
42 (declaim (inline two-sum
))
44 "Computes fl(a+b) and err(a+b)"
45 (declare (double-float a b
))
52 (declaim (inline two-prod
))
53 (declaim (inline split
))
54 ;; This algorithm is the version given by Yozo Hida. It has problems
55 ;; with overflow because we multiply by 1+2^27.
57 ;; But be very careful about replacing this with a new algorithm. The
58 ;; values computed here are very important to get the rounding right.
59 ;; If you change this, the rounding may be different, which will
60 ;; affect other parts of the algorithm.
62 ;; I (rtoy) tried a different algorithm that split the number in two
63 ;; as described, but without overflow. However, that caused
64 ;; -9.4294948327242751340284975915175w0/1w14 to return a value that
65 ;; wasn't really close to -9.4294948327242751340284975915175w-14.
67 ;; This also means we can't print numbers like 1w308 with the current
68 ;; printing algorithm, or even divide 1w308 by 10.
71 "Split the double-float number a into a-hi and a-lo such that a =
72 a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
73 a-lo contains the lower 26 bits."
74 (declare (double-float a
))
75 (let* ((tmp (* a
(+ 1 (expt 2 27))))
76 (a-hi (- tmp
(- tmp a
)))
82 "Split the double-float number a into a-hi and a-lo such that a =
83 a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
84 a-lo contains the lower 26 bits."
85 (declare (double-float a
)
87 ;; This splits the number a into 2 halves of 26 bits each, but the
88 ;; halves are, I think, supposed to be properly rounded in an IEEE
91 ;; For numbers that are very large, we use a different algorithm.
92 ;; For smaller numbers, we can use the original algorithm of Yozo
94 (if (> (abs a
) (scale-float 1d0
(- 1023 27)))
95 ;; I've tested this algorithm against Yozo's method for 1
96 ;; billion randomly generated double-floats between 2^(-995) and
97 ;; 2^996, and identical results are obtained. For numbers that
98 ;; are very small, this algorithm produces different numbers
99 ;; because of underflow. For very large numbers, we, of course
100 ;; produce different results because Yozo's method causes
102 (let* ((tmp (* a
(+ 1 (scale-float 1d0 -
27))))
103 (as (* a
(scale-float 1d0 -
27)))
104 (a-hi (* (- tmp
(- tmp as
)) (expt 2 27)))
108 (let* ((tmp (* a
(+ 1 (expt 2 27))))
109 (a-hi (- tmp
(- tmp a
)))
111 (values a-hi a-lo
))))
114 (defun two-prod (a b
)
115 "Compute fl(a*b) and err(a*b)"
116 (declare (double-float a b
))
118 (multiple-value-bind (a-hi a-lo
)
120 (multiple-value-bind (b-hi b-lo
)
122 (let ((e (+ (+ (- (* a-hi b-hi
) p
)
128 (declaim (inline two-sqr
))
130 "Compute fl(a*a) and err(a*b). This is a more efficient
131 implementation of two-prod"
132 (declare (double-float a
))
134 (multiple-value-bind (a-hi a-lo
)
136 (values q
(+ (+ (- (* a-hi a-hi
) q
)