qd-dd.lisp 8.03 KB
Newer Older
Raymond Toy's avatar
Raymond Toy committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
;;;; -*- Mode: lisp -*-
;;;;
;;;; Copyright (c) 2007 Raymond Toy
;;;;
;;;; Permission is hereby granted, free of charge, to any person
;;;; obtaining a copy of this software and associated documentation
;;;; files (the "Software"), to deal in the Software without
;;;; restriction, including without limitation the rights to use,
;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
;;;; copies of the Software, and to permit persons to whom the
;;;; Software is furnished to do so, subject to the following
;;;; conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
25

26
(in-package #:octi)
27
28
29
30

;;; double-double float routines needed for quad-double.
;;;
;;; Not needed for CMUCL.
31
32
;;;
;;; These routines were taken directly from CMUCL.
33

34
#||
35
36
37
(declaim (inline quick-two-sum))
(defun quick-two-sum (a b)
  "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
38
  (declare (double-float a b)
39
	   (optimize (speed 3) (safety 0) (debug 0)))
40
41
  (let* ((s (+ a b))
	 (e (- b (- s a))))
42
    (declare (double-float s e))
43
    (values s e)))
44
45
46
47
48
49
50
51
52
53
||#

(defmacro quick-two-sum (s e x y)
  (let ((a (gensym))
	(b (gensym)))
    `(let* ((,a ,x)
	    (,b ,y))
      (declare (double-float ,a ,b ,s ,e))
      (setf ,s (+ ,a ,b))
      (setf ,e (- ,b (- ,s ,a))))))
54

55
#||
56
57
58
(declaim (inline two-sum))
(defun two-sum (a b)
  "Computes fl(a+b) and err(a+b)"
59
  (declare (double-float a b)
60
	   (optimize (speed 3) (safety 0) (debug 0)))
61
62
63
64
  (let* ((s (+ a b))
	 (v (- s a))
	 (e (+ (- a (- s v))
	       (- b v))))
65
    (declare (double-float s v e))
66
    (values s e)))
67
||#
68

69
70
71
72
(defmacro two-sum (s e x y)
  "Computes fl(a+b) and err(a+b)"
  (let ((a (gensym))
	(b (gensym))
73
	(v (gensym)))
74
75
76
77
78
    `(let ((,a ,x)
	   (,b ,y))
      (declare (double-float ,a ,b))
      (setf ,s (+ ,a ,b))
      (let ((,v (- ,s ,a)))
79
80
81
	(declare (double-float ,v))
	(setf ,e (+ (- ,a (- ,s ,v))
		    (- ,b ,v)))))))
82

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
(declaim (inline two-prod))
(declaim (inline split))
;; This algorithm is the version given by Yozo Hida.  It has problems
;; with overflow because we multiply by 1+2^27.
;;
;; But be very careful about replacing this with a new algorithm.  The
;; values computed here are very important to get the rounding right.
;; If you change this, the rounding may be different, which will
;; affect other parts of the algorithm.
;;
;; I (rtoy) tried a different algorithm that split the number in two
;; as described, but without overflow.  However, that caused
;; -9.4294948327242751340284975915175w0/1w14 to return a value that
;; wasn't really close to -9.4294948327242751340284975915175w-14.
;;
;; This also means we can't print numbers like 1w308 with the current
;; printing algorithm, or even divide 1w308 by 10.
#+nil
(defun split (a)
  "Split the double-float number a into a-hi and a-lo such that a =
  a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
  a-lo contains the lower 26 bits."
  (declare (double-float a))
  (let* ((tmp (* a (+ 1 (expt 2 27))))
	 (a-hi (- tmp (- tmp a)))
	 (a-lo (- a a-hi)))
    (values a-hi a-lo)))


112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
;; Here is my limited understanding of what SPLIT is really supposed
;; to do.
;;
;; Let A be a double-float number.  We want to split the fraction bits
;; into 2 parts of 26 bits each.  This is best explained by example.
;; Let use use A = 1d50.  Use INTEGER-DECODE-FLOAT to display the bits
;; of A:
;;
;; (write (integer-decode-float 1d50) :base 2) ->
;;   10001000110110000111011000101011111100110010010011010
;;
;; Break this into 2 parts with the lower part having 27 bits and the
;; upper having 26 (because of the extra "hidden" bit):
;;
;;   10001000110110000111011000 101011111100110010010011010
;;
;; But this is not enough.  Note that the bottom half has a leading 1.
;; We want to round up the upper part.  Then we need to account for
;; this in the lower part:
;;
;;   10001000110110000111011001 -10100000011001101101100110
;;
;; This is the answer we want.  Convert these back to floats with the
;; appropriate exponents, and we get:
;;
;; 1.0000000087331024d50 and -8.733102285997912d41
;;
;; While this example worked out nicely, we should note that the
;; rounding operation above should be done in an IEEE round-to-even
;; fashion.  So if the lower part of the bits is exactly "half", we
;; round the upper part to even.  Thus,
;;
;; (float #b10001000110110000111011000100000000000000000000000000 1d0)
;; should be split into two parts:
;;
;;   10001000110110000111011000 100000000000000000000000000
;;
;; but
;;
;; (float #b10001000110110000111011001100000000000000000000000000 1d0)
;; is
;;
;;   10001000110110000111011010 -100000000000000000000000000

156
157
158
159
160
161
162
163
164
165
166
167
168
(defun split (a)
  "Split the double-float number a into a-hi and a-lo such that a =
  a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
  a-lo contains the lower 26 bits."
  (declare (double-float a)
	   (optimize (speed 3)))
  ;; This splits the number a into 2 halves of 26 bits each, but the
  ;; halves are, I think, supposed to be properly rounded in an IEEE
  ;; fashion.
  ;;
  ;; For numbers that are very large, we use a different algorithm.
  ;; For smaller numbers, we can use the original algorithm of Yozo
  ;; Hida.
169
  (if (> (abs a) #.(scale-float 1d0 (- 1023 27)))
170
171
172
173
174
175
176
      ;; I've tested this algorithm against Yozo's method for 1
      ;; billion randomly generated double-floats between 2^(-995) and
      ;; 2^996, and identical results are obtained.  For numbers that
      ;; are very small, this algorithm produces different numbers
      ;; because of underflow.  For very large numbers, we, of course
      ;; produce different results because Yozo's method causes
      ;; overflow.
177
178
179
      (let* ((tmp (* a #.(+ 1 (scale-float 1d0 -27))))
	     (as (* a #.(scale-float 1d0 -27)))
	     (a-hi (* (- tmp (- tmp as)) #.(scale-float 1d0 27)))
180
	     (a-lo (- a a-hi)))
181
	(declare (double-float tmp as a-hi a-lo))
182
183
	(values a-hi a-lo))
      ;; Yozo's algorithm.
184
      (let* ((tmp (* a #.(float (+ 1 (expt 2 27)) 1d0)))
185
186
	     (a-hi (- tmp (- tmp a)))
	     (a-lo (- a a-hi)))
187
	(declare (double-float tmp a-hi a-lo))
188
189
190
	(values a-hi a-lo))))


191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
;; Note that if you have an architecture that has a fused
;; multiply-subtract instruction that computes a*b-c with exactly one
;; rounding operation, you can use that instead of the complicated
;; routine below.  Power PC chips have such an instruction.
;;
;; Here is the code to do that, where (fused-multiply-subtract a b p)
;; computes a*b-p.
;;
;; (defun two-prod (a b)
;;   "Compute fl(a*b) and err(a*b)"
;;   (declare (double-float a b))
;;   (let* ((p (* a b))
;; 	 (err (fused-multiply-subtract a b p)))
;;     (values p err)))

206
207
(defun two-prod (a b)
  "Compute fl(a*b) and err(a*b)"
208
  (declare (double-float a b)
209
	   (optimize (speed 3) (safety #-allegro 0 #+allegro 1) (debug 0)))
210
  (let ((p (* a b)))
211
    (declare (double-float p))
212
213
    (multiple-value-bind (a-hi a-lo)
	(split a)
214
      (declare (double-float a-hi a-lo))
215
216
      (multiple-value-bind (b-hi b-lo)
	  (split b)
217
	(declare (double-float b-hi b-lo))
218
219
220
221
	(let ((e (+ (+ (- (* a-hi b-hi) p)
		       (* a-hi b-lo)
		       (* a-lo b-hi))
		    (* a-lo b-lo))))
222
	  (declare (double-float e))
223
224
225
226
227
228
	  (values p e))))))

(declaim (inline two-sqr))
(defun two-sqr (a)
  "Compute fl(a*a) and err(a*b).  This is a more efficient
  implementation of two-prod"
229
  (declare (double-float a)
230
	   (optimize (speed 3) (safety #-allegro 0 #+allegro 1) (debug 0)))
231
  (let ((q (* a a)))
232
    (declare (double-float q))
233
234
    (multiple-value-bind (a-hi a-lo)
	(split a)
235
      (declare (double-float a-hi a-lo))
236
237
238
      (values q (+ (+ (- (* a-hi a-hi) q)
		      (* 2 a-hi a-lo))
		   (* a-lo a-lo))))))