1 | ;;;; -*- Mode: lisp -*- |
---|
2 | ;;;; |
---|
3 | ;;;; Copyright (c) 2007, 2008, 2011 Raymond Toy |
---|
4 | ;;;; |
---|
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 |
---|
12 | ;;;; conditions: |
---|
13 | ;;;; |
---|
14 | ;;;; The above copyright notice and this permission notice shall be |
---|
15 | ;;;; included in all copies or substantial portions of the Software. |
---|
16 | ;;;; |
---|
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. |
---|
25 | |
---|
26 | ;;; This file contains the core routines for basic arithmetic |
---|
27 | ;;; operations on a %quad-double. This includes addition, |
---|
28 | ;;; subtraction, multiplication, division, negation. and absolute |
---|
29 | ;;; value. Some additional helper functions are included such as |
---|
30 | ;;; raising to an integer power. and the n'th root of a (non-negative) |
---|
31 | ;;; %quad-double. The basic comparison operators are included, and |
---|
32 | ;;; some simple tests for zerop, onep, plusp, and minusp. |
---|
33 | ;;; |
---|
34 | ;;; The basic algorithms are based on Yozo Hida's double-double |
---|
35 | ;;; implementation. However, some were copied from CMUCL and modified |
---|
36 | ;;; to support quad-doubles. |
---|
37 | |
---|
38 | |
---|
39 | (in-package #:octi) |
---|
40 | |
---|
41 | #+cmu |
---|
42 | (eval-when (:compile-toplevel) |
---|
43 | (setf ext:*inline-expansion-limit* 1600)) |
---|
44 | |
---|
45 | ;; All of the following functions should be inline. |
---|
46 | ;;(declaim (inline three-sum2)) |
---|
47 | |
---|
48 | ;; Internal routines for implementing quad-double. |
---|
49 | |
---|
50 | #+cmu |
---|
51 | (defmacro two-sum (s e x y) |
---|
52 | `(multiple-value-setq (,s ,e) |
---|
53 | (c::two-sum ,x ,y))) |
---|
54 | |
---|
55 | |
---|
56 | (defmacro three-sum (s0 s1 s2 a b c) |
---|
57 | (let ((aa (gensym)) |
---|
58 | (bb (gensym)) |
---|
59 | (cc (gensym)) |
---|
60 | (t1 (gensym)) |
---|
61 | (t2 (gensym)) |
---|
62 | (t3 (gensym))) |
---|
63 | `(let* ((,aa ,a) |
---|
64 | (,bb ,b) |
---|
65 | (,cc ,c) |
---|
66 | (,t1 0d0) |
---|
67 | (,t2 ,t1) |
---|
68 | (,t3 ,t1)) |
---|
69 | (declare (double-float ,t1 ,t2 ,t3)) |
---|
70 | (two-sum ,t1 ,t2 ,aa ,bb) |
---|
71 | (two-sum ,s0 ,t3 ,cc ,t1) |
---|
72 | (two-sum ,s1 ,s2 ,t2 ,t3)))) |
---|
73 | |
---|
74 | |
---|
75 | |
---|
76 | #+nil |
---|
77 | (defun three-sum2 (a b c) |
---|
78 | (declare (double-float a b c) |
---|
79 | (optimize (speed 3))) |
---|
80 | (let* ((t1 0d0) |
---|
81 | (t2 t1) |
---|
82 | (t3 t1)) |
---|
83 | (two-sum t1 t2 a b) |
---|
84 | (two-sum a t3 c t1) |
---|
85 | (values a (cl:+ t2 t3)))) |
---|
86 | |
---|
87 | (defmacro three-sum2 (s0 s1 a b c) |
---|
88 | (let ((aa (gensym)) |
---|
89 | (bb (gensym)) |
---|
90 | (cc (gensym)) |
---|
91 | (t1 (gensym)) |
---|
92 | (t2 (gensym)) |
---|
93 | (t3 (gensym))) |
---|
94 | `(let* ((,aa ,a) |
---|
95 | (,bb ,b) |
---|
96 | (,cc ,c) |
---|
97 | (,t1 0d0) |
---|
98 | (,t2 ,t1) |
---|
99 | (,t3 ,t1)) |
---|
100 | (declare (double-float ,t1 ,t2 ,t3)) |
---|
101 | (two-sum ,t1 ,t2 ,aa ,bb) |
---|
102 | (two-sum ,s0 ,t3 ,cc ,t1) |
---|
103 | (setf ,s1 (+ ,t2 ,t3))))) |
---|
104 | |
---|
105 | ;; Not needed???? |
---|
106 | #+nil |
---|
107 | (declaim (inline quick-three-accum)) |
---|
108 | #+nil |
---|
109 | (defun quick-three-accum (a b c) |
---|
110 | (declare (double-float a b c) |
---|
111 | (optimize (speed 3) (space 0))) |
---|
112 | (multiple-value-bind (s b) |
---|
113 | (two-sum b c) |
---|
114 | (multiple-value-bind (s a) |
---|
115 | (two-sum a s) |
---|
116 | (let ((za (/= a 0)) |
---|
117 | (zb (/= b 0))) |
---|
118 | (when (and za zb) |
---|
119 | (return-from quick-three-accum (values (cl:+ s 0d0) (cl:+ a 0d0) (cl:+ b 0d0)))) |
---|
120 | (when (not za) |
---|
121 | (setf a s) |
---|
122 | (setf s 0d0)) |
---|
123 | (when (not zb) |
---|
124 | (setf b a) |
---|
125 | (setf a s)) |
---|
126 | (values 0d0 a b))))) |
---|
127 | |
---|
128 | |
---|
129 | ;; These functions are quite short, so we inline them to minimize |
---|
130 | ;; consing. |
---|
131 | (declaim (inline make-qd-d |
---|
132 | add-d-qd |
---|
133 | add-dd-qd |
---|
134 | neg-qd |
---|
135 | neg-qd-t |
---|
136 | sub-qd |
---|
137 | sub-qd-dd |
---|
138 | sub-qd-d |
---|
139 | sub-d-qd |
---|
140 | #+cmu |
---|
141 | make-qd-dd |
---|
142 | abs-qd |
---|
143 | qd-> |
---|
144 | qd-< |
---|
145 | qd->= |
---|
146 | qd-<= |
---|
147 | zerop-qd |
---|
148 | onep-qd |
---|
149 | plusp-qd |
---|
150 | minusp-qd |
---|
151 | qd-= |
---|
152 | scale-float-qd)) |
---|
153 | |
---|
154 | ;; Should these functions be inline? The QD C++ source has these as |
---|
155 | ;; inline functions, but these are fairly large functions. However, |
---|
156 | ;; inlining makes quite a big difference in speed and consing. |
---|
157 | #+cmu |
---|
158 | (declaim (#+qd-inline inline |
---|
159 | #-qd-inline ext:maybe-inline |
---|
160 | renorm-4 |
---|
161 | renorm-5 |
---|
162 | add-qd-d |
---|
163 | add-qd-d-t |
---|
164 | add-qd-dd |
---|
165 | add-qd |
---|
166 | add-qd-t |
---|
167 | mul-qd-d |
---|
168 | mul-qd-d-t |
---|
169 | mul-qd-dd |
---|
170 | mul-qd |
---|
171 | mul-qd-t |
---|
172 | sqr-qd |
---|
173 | sqr-qd-t |
---|
174 | div-qd |
---|
175 | div-qd-t |
---|
176 | div-qd-d |
---|
177 | div-qd-dd)) |
---|
178 | |
---|
179 | #+cmu |
---|
180 | (defmacro quick-two-sum (s e x y) |
---|
181 | `(multiple-value-setq (,s ,e) |
---|
182 | (c::quick-two-sum ,x ,y))) |
---|
183 | |
---|
184 | #-(or qd-inline (not cmu)) |
---|
185 | (declaim (ext:start-block renorm-4 renorm-5 |
---|
186 | make-qd-d |
---|
187 | add-qd-d add-d-qd add-qd-dd |
---|
188 | add-dd-qd |
---|
189 | add-qd add-qd-t |
---|
190 | add-qd-d-t |
---|
191 | neg-qd |
---|
192 | sub-qd sub-qd-dd sub-qd-d sub-d-qd |
---|
193 | mul-qd-d mul-qd-dd |
---|
194 | mul-qd |
---|
195 | mul-qd-t |
---|
196 | mul-qd-d-t |
---|
197 | sqr-qd |
---|
198 | sqr-qd-t |
---|
199 | div-qd div-qd-d div-qd-dd |
---|
200 | div-qd-t |
---|
201 | make-qd-dd |
---|
202 | )) |
---|
203 | |
---|
204 | #+(or) |
---|
205 | (defun quick-renorm (c0 c1 c2 c3 c4) |
---|
206 | (declare (double-float c0 c1 c2 c3 c4) |
---|
207 | (optimize (speed 3))) |
---|
208 | (multiple-value-bind (s t3) |
---|
209 | (quick-two-sum c3 c4) |
---|
210 | (multiple-value-bind (s t2) |
---|
211 | (quick-two-sum c2 s) |
---|
212 | (multiple-value-bind (s t1) |
---|
213 | (quick-two-sum c1 s) |
---|
214 | (multiple-value-bind (c0 t0) |
---|
215 | (quick-two-sum c0 s) |
---|
216 | (multiple-value-bind (s t2) |
---|
217 | (quick-two-sum t2 t3) |
---|
218 | (multiple-value-bind (s t1) |
---|
219 | (quick-two-sum t1 s) |
---|
220 | (multiple-value-bind (c1 t0) |
---|
221 | (quick-two-sum t0 s) |
---|
222 | (multiple-value-bind (s t1) |
---|
223 | (quick-two-sum t1 t2) |
---|
224 | (multiple-value-bind (c2 t0) |
---|
225 | (quick-two-sum t0 s) |
---|
226 | (values c0 c1 c2 (cl:+ t0 t1)))))))))))) |
---|
227 | |
---|
228 | (defun renorm-4 (c0 c1 c2 c3) |
---|
229 | (declare (double-float c0 c1 c2 c3) |
---|
230 | (optimize (speed 3) (safety #-allegro 0 #+allegro 1) (debug 0))) |
---|
231 | (let ((s0 0d0) |
---|
232 | (s1 0d0) |
---|
233 | (s2 0d0) |
---|
234 | (s3 0d0)) |
---|
235 | (declare (double-float s0 s1 s2 s3)) |
---|
236 | (quick-two-sum s0 c3 c2 c3) |
---|
237 | (quick-two-sum s0 c2 c1 s0) |
---|
238 | (quick-two-sum c0 c1 c0 s0) |
---|
239 | (setf s0 c0) |
---|
240 | (setf s1 c1) |
---|
241 | (cond ((/= s1 0) |
---|
242 | (quick-two-sum s1 s2 s1 c2) |
---|
243 | (if (/= s2 0) |
---|
244 | (quick-two-sum s2 s3 s2 c3) |
---|
245 | (quick-two-sum s1 s2 s1 c3))) |
---|
246 | (t |
---|
247 | (quick-two-sum s0 s1 s0 c2) |
---|
248 | (if (/= s1 0) |
---|
249 | (quick-two-sum s1 s2 s1 c3) |
---|
250 | (quick-two-sum s0 s1 s0 c3)))) |
---|
251 | (values s0 s1 s2 s3))) |
---|
252 | |
---|
253 | (defun renorm-5 (c0 c1 c2 c3 c4) |
---|
254 | (declare (double-float c0 c1 c2 c3 c4) |
---|
255 | (optimize (speed 3) (safety #-allegro 0 #+allegro 1))) |
---|
256 | (let ((s0 0d0) |
---|
257 | (s1 0d0) |
---|
258 | (s2 0d0) |
---|
259 | (s3 0d0)) |
---|
260 | (declare (double-float s0 s1 s2 s3)) |
---|
261 | (quick-two-sum s0 c4 c3 c4) |
---|
262 | (quick-two-sum s0 c3 c2 s0) |
---|
263 | (quick-two-sum s0 c2 c1 s0) |
---|
264 | (quick-two-sum c0 c1 c0 s0) |
---|
265 | (quick-two-sum s0 s1 c0 c1) |
---|
266 | (cond ((/= s1 0) |
---|
267 | (quick-two-sum s1 s2 s1 c2) |
---|
268 | (cond ((/= s2 0) |
---|
269 | (quick-two-sum s2 s3 s2 c3) |
---|
270 | (if (/= s3 0) |
---|
271 | (incf s3 c4) |
---|
272 | (incf s2 c4))) |
---|
273 | (t |
---|
274 | (quick-two-sum s1 s2 s1 c3) |
---|
275 | (if (/= s2 0) |
---|
276 | (quick-two-sum s2 s3 s2 c4) |
---|
277 | (quick-two-sum s1 s2 s1 c4))))) |
---|
278 | (t |
---|
279 | (quick-two-sum s0 s1 s0 c2) |
---|
280 | (cond ((/= s1 0) |
---|
281 | (quick-two-sum s1 s2 s1 c3) |
---|
282 | (if (/= s2 0) |
---|
283 | (quick-two-sum s2 s3 s2 c4) |
---|
284 | (quick-two-sum s1 s2 s1 c4))) |
---|
285 | (t |
---|
286 | (quick-two-sum s0 s1 s0 c3) |
---|
287 | (if (/= s1 0) |
---|
288 | (quick-two-sum s1 s2 s1 c4) |
---|
289 | (quick-two-sum s0 s1 s0 c4)))))) |
---|
290 | (values s0 s1 s2 s3))) |
---|
291 | |
---|
292 | (defun make-qd-d (a0 &optional (a1 0d0 a1-p) (a2 0d0) (a3 0d0)) |
---|
293 | "Create a %QUAD-DOUBLE from four double-floats, appropriately |
---|
294 | normalizing the result from the four double-floats. A0 is the most |
---|
295 | significant part of the %QUAD-DOUBLE, and A3 is the least. The optional |
---|
296 | parameters default to 0. |
---|
297 | " |
---|
298 | (declare (double-float a0 a1 a2 a3) |
---|
299 | (optimize (speed 3) |
---|
300 | #+cmu |
---|
301 | (ext:inhibit-warnings 3))) |
---|
302 | (if a1-p |
---|
303 | (multiple-value-bind (s0 s1 s2 s3) |
---|
304 | (renorm-4 a0 a1 a2 a3) |
---|
305 | (%make-qd-d s0 s1 s2 s3)) |
---|
306 | (%make-qd-d a0 0d0 0d0 0d0))) |
---|
307 | |
---|
308 | ;;;; Addition |
---|
309 | |
---|
310 | ;; Quad-double + double |
---|
311 | (defun add-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
312 | "Return the sum of the %QUAD-DOUBLE A and the DOUBLE-FLOAT B. |
---|
313 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
314 | (add-qd-d-t a b target)) |
---|
315 | |
---|
316 | (defun add-qd-d-t (a b target) |
---|
317 | "Add a quad-double A and a double-float B" |
---|
318 | (declare (type %quad-double a #+oct-array target) |
---|
319 | (double-float b) |
---|
320 | (optimize (speed 3) |
---|
321 | (space 0)) |
---|
322 | (inline float-infinity-p) |
---|
323 | #+(and cmu (not oct-array)) (ignore target)) |
---|
324 | (let* ((c0 0d0) |
---|
325 | (e c0) |
---|
326 | (c1 c0) |
---|
327 | (c2 c0) |
---|
328 | (c3 c0)) |
---|
329 | (declare (double-float e c0 c1 c2 c3)) |
---|
330 | (two-sum c0 e (qd-0 a) b) |
---|
331 | |
---|
332 | (when (float-infinity-p c0) |
---|
333 | (return-from add-qd-d-t (%store-qd-d target c0 0d0 0d0 0d0))) |
---|
334 | (two-sum c1 e (qd-1 a) e) |
---|
335 | (two-sum c2 e (qd-2 a) e) |
---|
336 | (two-sum c3 e (qd-3 a) e) |
---|
337 | (multiple-value-bind (r0 r1 r2 r3) |
---|
338 | (renorm-5 c0 c1 c2 c3 e) |
---|
339 | (if (and (zerop (qd-0 a)) (zerop b)) |
---|
340 | (%store-qd-d target c0 0d0 0d0 0d0) |
---|
341 | (%store-qd-d target r0 r1 r2 r3))))) |
---|
342 | |
---|
343 | (defun add-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
344 | "Return the sum of the DOUBLE-FLOAT A and the %QUAD-DOUBLE B. |
---|
345 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
346 | (declare (double-float a) |
---|
347 | (type %quad-double b) |
---|
348 | (optimize (speed 3)) |
---|
349 | #+(and cmu (not oct-array)) (ignore target)) |
---|
350 | (add-qd-d b a #+oct-array target)) |
---|
351 | |
---|
352 | #+cmu |
---|
353 | (defun add-qd-dd (a b) |
---|
354 | "Add a quad-double A and a double-double B" |
---|
355 | (declare (type %quad-double a) |
---|
356 | (double-double-float b) |
---|
357 | (optimize (speed 3) |
---|
358 | (space 0))) |
---|
359 | (let* ((s0 0d0) |
---|
360 | (t0 s0) |
---|
361 | (s1 s0) |
---|
362 | (t1 s0) |
---|
363 | (s2 s0) |
---|
364 | (s3 s0)) |
---|
365 | (declare (double-float s0 s1 s3 t0 t1 s2)) |
---|
366 | (two-sum s0 t1 (qd-0 a) (kernel:double-double-hi b)) |
---|
367 | (two-sum s1 t1 (qd-1 a) (kernel:double-double-lo b)) |
---|
368 | (two-sum s1 t0 s1 t0) |
---|
369 | (three-sum s2 t0 t1 (qd-2 a) t0 t1) |
---|
370 | (two-sum s3 t0 t0 (qd-3 a)) |
---|
371 | (incf t0 t1) |
---|
372 | (multiple-value-bind (a0 a1 a2 a3) |
---|
373 | (renorm-5 s0 s1 s2 s3 t0) |
---|
374 | (declare (double-float a0 a1 a2 a3)) |
---|
375 | (%make-qd-d a0 a1 a2 a3)))) |
---|
376 | |
---|
377 | #+cmu |
---|
378 | (defun add-dd-qd (a b) |
---|
379 | (declare (double-double-float a) |
---|
380 | (type %quad-double b) |
---|
381 | (optimize (speed 3) |
---|
382 | (space 0))) |
---|
383 | (add-qd-dd b a)) |
---|
384 | |
---|
385 | |
---|
386 | #+(or) |
---|
387 | (defun add-qd-1 (a b) |
---|
388 | (declare (type %quad-double a b) |
---|
389 | (optimize (speed 3))) |
---|
390 | (multiple-value-bind (s0 t0) |
---|
391 | (two-sum (qd-0 a) (qd-0 b)) |
---|
392 | (multiple-value-bind (s1 t1) |
---|
393 | (two-sum (qd-1 a) (qd-1 b)) |
---|
394 | (multiple-value-bind (s2 t2) |
---|
395 | (two-sum (qd-2 a) (qd-2 b)) |
---|
396 | (multiple-value-bind (s3 t3) |
---|
397 | (two-sum (qd-3 a) (qd-3 b)) |
---|
398 | (multiple-value-bind (s1 t0) |
---|
399 | (two-sum s1 t0) |
---|
400 | (multiple-value-bind (s2 t0 t1) |
---|
401 | (three-sum s2 t0 t1) |
---|
402 | (multiple-value-bind (s3 t0) |
---|
403 | (three-sum2 s3 t0 t2) |
---|
404 | (let ((t0 (cl:+ t0 t1 t3))) |
---|
405 | (multiple-value-call #'%make-qd-d |
---|
406 | (renorm-5 s0 s1 s2 s3 t0))))))))))) |
---|
407 | |
---|
408 | ;; Same as above, except that everything is expanded out for compilers |
---|
409 | ;; which don't do a very good job with dataflow. CMUCL is one of |
---|
410 | ;; those compilers. |
---|
411 | |
---|
412 | (defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
413 | "Return the sum of the %QUAD-DOUBLE numbers A and B. |
---|
414 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
415 | (add-qd-t a b target)) |
---|
416 | |
---|
417 | |
---|
418 | (defun add-qd-t (a b target) |
---|
419 | (declare (type %quad-double a b #+oct-array target) |
---|
420 | (optimize (speed 3) |
---|
421 | (space 0)) |
---|
422 | #+(and cmu (not oct-array)) |
---|
423 | (ignore target)) |
---|
424 | ;; This is the version that is NOT IEEE. Should we use the IEEE |
---|
425 | ;; version? It's quite a bit more complicated. |
---|
426 | ;; |
---|
427 | ;; In addition, this is reorganized to minimize data dependency. |
---|
428 | (with-qd-parts (a0 a1 a2 a3) |
---|
429 | a |
---|
430 | (declare (double-float a0 a1 a2 a3)) |
---|
431 | (with-qd-parts (b0 b1 b2 b3) |
---|
432 | b |
---|
433 | (declare (double-float b0 b1 b2 b3)) |
---|
434 | (let ((s0 (cl:+ a0 b0)) |
---|
435 | (s1 (cl:+ a1 b1)) |
---|
436 | (s2 (cl:+ a2 b2)) |
---|
437 | (s3 (cl:+ a3 b3))) |
---|
438 | (declare (double-float s0 s1 s2 s3) |
---|
439 | (inline float-infinity-p)) |
---|
440 | |
---|
441 | (when (float-infinity-p s0) |
---|
442 | (return-from add-qd-t (%store-qd-d target s0 0d0 0d0 0d0))) |
---|
443 | (let ((v0 (cl:- s0 a0)) |
---|
444 | (v1 (cl:- s1 a1)) |
---|
445 | (v2 (cl:- s2 a2)) |
---|
446 | (v3 (cl:- s3 a3))) |
---|
447 | (declare (double-float v0 v1 v2 v3)) |
---|
448 | (let ((u0 (cl:- s0 v0)) |
---|
449 | (u1 (cl:- s1 v1)) |
---|
450 | (u2 (cl:- s2 v2)) |
---|
451 | (u3 (cl:- s3 v3))) |
---|
452 | (declare (double-float u0 u1 u2 u3)) |
---|
453 | (let ((w0 (cl:- a0 u0)) |
---|
454 | (w1 (cl:- a1 u1)) |
---|
455 | (w2 (cl:- a2 u2)) |
---|
456 | (w3 (cl:- a3 u3))) |
---|
457 | (declare (double-float w0 w1 w2 w3)) |
---|
458 | (let ((u0 (cl:- b0 v0)) |
---|
459 | (u1 (cl:- b1 v1)) |
---|
460 | (u2 (cl:- b2 v2)) |
---|
461 | (u3 (cl:- b3 v3))) |
---|
462 | (declare (double-float u0 u1 u2 u3)) |
---|
463 | (let ((t0 (cl:+ w0 u0)) |
---|
464 | (t1 (cl:+ w1 u1)) |
---|
465 | (t2 (cl:+ w2 u2)) |
---|
466 | (t3 (cl:+ w3 u3))) |
---|
467 | (declare (double-float t0 t1 t2 t3)) |
---|
468 | (two-sum s1 t0 s1 t0) |
---|
469 | (three-sum s2 t0 t1 s2 t0 t1) |
---|
470 | (three-sum2 s3 t0 s3 t0 t2) |
---|
471 | (setf t0 (cl:+ t0 t1 t3)) |
---|
472 | ;; Renormalize |
---|
473 | (multiple-value-setq (s0 s1 s2 s3) |
---|
474 | (renorm-5 s0 s1 s2 s3 t0)) |
---|
475 | (if (and (zerop a0) (zerop b0)) |
---|
476 | (%store-qd-d target (+ a0 b0) 0d0 0d0 0d0) |
---|
477 | (%store-qd-d target s0 s1 s2 s3))))))))))) |
---|
478 | |
---|
479 | (defun neg-qd-t (a target) |
---|
480 | (declare (type %quad-double a #+oct-array target) |
---|
481 | #+(and cmu (not oct-array)) (ignore target)) |
---|
482 | (with-qd-parts (a0 a1 a2 a3) |
---|
483 | a |
---|
484 | (declare (double-float a0 a1 a2 a3)) |
---|
485 | (%store-qd-d target (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3)))) |
---|
486 | |
---|
487 | (defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
488 | "Return the negative of the %QUAD-DOUBLE number A. |
---|
489 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
490 | (neg-qd-t a target)) |
---|
491 | |
---|
492 | |
---|
493 | |
---|
494 | (defun sub-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
495 | "Return the difference between the %QUAD-DOUBLE numbers A and B. |
---|
496 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
497 | (declare (type %quad-double a b)) |
---|
498 | (add-qd-t a (neg-qd b) target)) |
---|
499 | |
---|
500 | (defun sub-qd-t (a b target) |
---|
501 | (add-qd-t a (neg-qd b) target)) |
---|
502 | |
---|
503 | #+cmu |
---|
504 | (defun sub-qd-dd (a b) |
---|
505 | (declare (type %quad-double a) |
---|
506 | (type double-double-float b)) |
---|
507 | (add-qd-dd a (cl:- b))) |
---|
508 | |
---|
509 | (defun sub-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
510 | "Return the difference between the %QUAD-DOUBLE number A and |
---|
511 | the DOUBLE-FLOAT number B. |
---|
512 | |
---|
513 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
514 | (declare (type %quad-double a) |
---|
515 | (type double-float b) |
---|
516 | #+(and cmu (not oct-array)) (ignore target)) |
---|
517 | (add-qd-d a (cl:- b) #+oct-array target)) |
---|
518 | |
---|
519 | (defun sub-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
520 | "Return the difference between the DOUBLE-FLOAT number A and |
---|
521 | the %QUAD-DOUBLE number B. |
---|
522 | |
---|
523 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
524 | (declare (type double-float a) |
---|
525 | (type %quad-double b) |
---|
526 | #+(and cmu (not oct-array)) (ignore target)) |
---|
527 | ;; a - b = a + (-b) |
---|
528 | (add-d-qd a (neg-qd b) #+oct-array target)) |
---|
529 | |
---|
530 | |
---|
531 | ;; Works |
---|
532 | ;; (mul-qd-d (sqrt-qd (make-qd-dd 2w0 0w0)) 10d0) -> |
---|
533 | ;; 14.1421356237309504880168872420969807856967187537694807317667973799q0 |
---|
534 | ;; |
---|
535 | ;; Clisp says |
---|
536 | ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 |
---|
537 | ;; |
---|
538 | (defun mul-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
539 | "Return the product of the %QUAD-DOUBLE number A and |
---|
540 | the DOUBLE-FLOAT number B. |
---|
541 | |
---|
542 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
543 | (mul-qd-d-t a b target)) |
---|
544 | |
---|
545 | (defun mul-qd-d-t (a b target) |
---|
546 | "Multiply quad-double A with B" |
---|
547 | (declare (type %quad-double a #+oct-array target) |
---|
548 | (double-float b) |
---|
549 | (optimize (speed 3) |
---|
550 | (space 0)) |
---|
551 | (inline float-infinity-p) |
---|
552 | #+(and cmu (not oct-array)) (ignore target)) |
---|
553 | (multiple-value-bind (p0 q0) |
---|
554 | (two-prod (qd-0 a) b) |
---|
555 | |
---|
556 | (when (float-infinity-p p0) |
---|
557 | (return-from mul-qd-d-t (%store-qd-d target p0 0d0 0d0 0d0))) |
---|
558 | (multiple-value-bind (p1 q1) |
---|
559 | (two-prod (qd-1 a) b) |
---|
560 | (declare (double-float p1 q1)) |
---|
561 | (multiple-value-bind (p2 q2) |
---|
562 | (two-prod (qd-2 a) b) |
---|
563 | (declare (double-float p2 q2)) |
---|
564 | (let* ((p3 (cl:* (qd-3 a) b)) |
---|
565 | (s0 p0) |
---|
566 | (s1 p0) |
---|
567 | (s2 p0)) |
---|
568 | (declare (double-float s0 s1 s2 p3)) |
---|
569 | (two-sum s1 s2 q0 p1) |
---|
570 | (three-sum s2 q1 p2 s2 q1 p2) |
---|
571 | (three-sum2 q1 q2 q1 q2 p3) |
---|
572 | (let ((s3 q1) |
---|
573 | (s4 (cl:+ q2 p2))) |
---|
574 | (multiple-value-bind (s0 s1 s2 s3) |
---|
575 | (renorm-5 s0 s1 s2 s3 s4) |
---|
576 | (if (zerop s0) |
---|
577 | (%store-qd-d target (float-sign p0 0d0) 0d0 0d0 0d0) |
---|
578 | (%store-qd-d target s0 s1 s2 s3))))))))) |
---|
579 | |
---|
580 | ;; a0 * b0 0 |
---|
581 | ;; a0 * b1 1 |
---|
582 | ;; a1 * b0 2 |
---|
583 | ;; a1 * b1 3 |
---|
584 | ;; a2 * b0 4 |
---|
585 | ;; a2 * b1 5 |
---|
586 | ;; a3 * b0 6 |
---|
587 | ;; a3 * b1 7 |
---|
588 | |
---|
589 | ;; Not working. |
---|
590 | ;; (mul-qd-dd (sqrt-qd (make-qd-dd 2w0 0w0)) 10w0) -> |
---|
591 | ;; 14.142135623730950488016887242097022172449805747901877456053837224q0 |
---|
592 | ;; |
---|
593 | ;; But clisp says |
---|
594 | ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 |
---|
595 | ;; ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
---|
596 | ;; |
---|
597 | ;; Running a test program using qd (2.1.210) shows that we get the |
---|
598 | ;; same wrong answer. |
---|
599 | #+(or) |
---|
600 | (defun mul-qd-dd (a b) |
---|
601 | (declare (type %quad-double a) |
---|
602 | (double-double-float b) |
---|
603 | (optimize (speed 3))) |
---|
604 | (multiple-value-bind (p0 q0) |
---|
605 | (two-prod (qd-0 a) (kernel:double-double-hi b)) |
---|
606 | (multiple-value-bind (p1 q1) |
---|
607 | (two-prod (qd-0 a) (kernel:double-double-lo b)) |
---|
608 | (multiple-value-bind (p2 q2) |
---|
609 | (two-prod (qd-1 a) (kernel:double-double-hi b)) |
---|
610 | (multiple-value-bind (p3 q3) |
---|
611 | (two-prod (qd-1 a) (kernel:double-double-lo b)) |
---|
612 | (multiple-value-bind (p4 q4) |
---|
613 | (two-prod (qd-2 a) (kernel:double-double-hi b)) |
---|
614 | (format t "p0, q0 = ~A ~A~%" p0 q0) |
---|
615 | (format t "p1, q1 = ~A ~A~%" p1 q1) |
---|
616 | (format t "p2, q2 = ~A ~A~%" p2 q2) |
---|
617 | (format t "p3, q3 = ~A ~A~%" p3 q3) |
---|
618 | (format t "p4, q4 = ~A ~A~%" p4 q4) |
---|
619 | (multiple-value-setq (p1 p2 q0) |
---|
620 | (three-sum p1 p2 q0)) |
---|
621 | (format t "p1 = ~A~%" p1) |
---|
622 | (format t "p2 = ~A~%" p2) |
---|
623 | (format t "q0 = ~A~%" q0) |
---|
624 | ;; five-three-sum |
---|
625 | (multiple-value-setq (p2 p3 p4) |
---|
626 | (three-sum p2 p3 p4)) |
---|
627 | (format t "p2 = ~A~%" p2) |
---|
628 | (format t "p3 = ~A~%" p3) |
---|
629 | (format t "p4 = ~A~%" p4) |
---|
630 | (multiple-value-setq (q1 q2) |
---|
631 | (two-sum q1 q2)) |
---|
632 | (multiple-value-bind (s0 t0) |
---|
633 | (two-sum p2 q1) |
---|
634 | (multiple-value-bind (s1 t1) |
---|
635 | (two-sum p3 q2) |
---|
636 | (multiple-value-setq (s1 t0) |
---|
637 | (two-sum s1 t0)) |
---|
638 | (let ((s2 (cl:+ t0 t1 p4)) |
---|
639 | (p2 s0) |
---|
640 | (p3 (cl:+ (cl:* (qd-2 a) |
---|
641 | (kernel:double-double-hi b)) |
---|
642 | (cl:* (qd-3 a) |
---|
643 | (kernel:double-double-lo b)) |
---|
644 | q3 q4))) |
---|
645 | (multiple-value-setq (p3 q0) |
---|
646 | (three-sum2 p3 q0 s1)) |
---|
647 | (let ((p4 (cl:+ q0 s2))) |
---|
648 | (multiple-value-call #'%make-qd-d |
---|
649 | (renorm-5 p0 p1 p2 p3 p4)))))))))))) |
---|
650 | |
---|
651 | ;; a0 * b0 0 |
---|
652 | ;; a0 * b1 1 |
---|
653 | ;; a1 * b0 2 |
---|
654 | ;; a0 * b2 3 |
---|
655 | ;; a1 * b1 4 |
---|
656 | ;; a2 * b0 5 |
---|
657 | ;; a0 * b3 6 |
---|
658 | ;; a1 * b2 7 |
---|
659 | ;; a2 * b1 8 |
---|
660 | ;; a3 * b0 9 |
---|
661 | |
---|
662 | ;; Works |
---|
663 | ;; (mul-qd (sqrt-qd (make-qd-dd 2w0 0w0)) (make-qd-dd 10w0 0w0)) -> |
---|
664 | ;; 14.1421356237309504880168872420969807856967187537694807317667973799q0 |
---|
665 | ;; |
---|
666 | ;; Clisp says |
---|
667 | ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 |
---|
668 | |
---|
669 | (defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
670 | "Returns the product of the %QUAD-DOUBLE numbers A and B. |
---|
671 | If TARGET is given, TARGET is destructively modified to contain the result." |
---|
672 | (mul-qd-t a b target)) |
---|
673 | |
---|
674 | (defun mul-qd-t (a b target) |
---|
675 | (declare (type %quad-double a b #+oct-array target) |
---|
676 | (optimize (speed 3) |
---|
677 | (space 0)) |
---|
678 | (inline float-infinity-p) |
---|
679 | #+(and cmu (not oct-array)) |
---|
680 | (ignore target)) |
---|
681 | (with-qd-parts (a0 a1 a2 a3) |
---|
682 | a |
---|
683 | (declare (double-float a0 a1 a2 a3)) |
---|
684 | (with-qd-parts (b0 b1 b2 b3) |
---|
685 | b |
---|
686 | (declare (double-float b0 b1 b2 b3)) |
---|
687 | (multiple-value-bind (p0 q0) |
---|
688 | (two-prod a0 b0) |
---|
689 | #+cmu |
---|
690 | (when (float-infinity-p p0) |
---|
691 | (return-from mul-qd-t (%store-qd-d target p0 0d0 0d0 0d0))) |
---|
692 | (multiple-value-bind (p1 q1) |
---|
693 | (two-prod a0 b1) |
---|
694 | (multiple-value-bind (p2 q2) |
---|
695 | (two-prod a1 b0) |
---|
696 | (multiple-value-bind (p3 q3) |
---|
697 | (two-prod a0 b2) |
---|
698 | (multiple-value-bind (p4 q4) |
---|
699 | (two-prod a1 b1) |
---|
700 | (multiple-value-bind (p5 q5) |
---|
701 | (two-prod a2 b0) |
---|
702 | ;; Start accumulation |
---|
703 | (three-sum p1 p2 q0 p1 p2 q0) |
---|
704 | |
---|
705 | ;; six-three-sum of p2, q1, q2, p3, p4, p5 |
---|
706 | (three-sum p2 q1 q2 p2 q1 q2) |
---|
707 | (three-sum p3 p4 p5 p3 p4 p5) |
---|
708 | ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5) |
---|
709 | (let* ((s0 0d0) |
---|
710 | (s1 s0) |
---|
711 | (t0 s0) |
---|
712 | (t1 s0)) |
---|
713 | (declare (double-float s0 s1 t0 t1)) |
---|
714 | (two-sum s0 t0 p2 p3) |
---|
715 | (two-sum s1 t1 q1 p4) |
---|
716 | (let ((s2 (cl:+ q2 p5))) |
---|
717 | (declare (double-float s2)) |
---|
718 | (two-sum s1 t0 s1 t0) |
---|
719 | (incf s2 (cl:+ t0 t1)) |
---|
720 | ;; O(eps^3) order terms. This is the sloppy |
---|
721 | ;; multiplication version. Should we use |
---|
722 | ;; the precise version? It's significantly |
---|
723 | ;; more complex. |
---|
724 | |
---|
725 | (incf s1 (cl:+ (cl:* a0 b3) |
---|
726 | (cl:* a1 b2) |
---|
727 | (cl:* a2 b1) |
---|
728 | (cl:* a3 b0) |
---|
729 | q0 q3 q4 q5)) |
---|
730 | #+nil |
---|
731 | (format t "p0,p1,s0,s1,s2 = ~a ~a ~a ~a ~a~%" |
---|
732 | p0 p1 s0 s1 s2) |
---|
733 | (multiple-value-bind (r0 r1 s0 s1) |
---|
734 | (renorm-5 p0 p1 s0 s1 s2) |
---|
735 | (if (zerop r0) |
---|
736 | (%store-qd-d target p0 0d0 0d0 0d0) |
---|
737 | (%store-qd-d target r0 r1 s0 s1)))))))))))))) |
---|
738 | |
---|
739 | |
---|
740 | ;; This is the non-sloppy version. I think this works just fine, but |
---|
741 | ;; since qd defaults to the sloppy multiplication version, we do the |
---|
742 | ;; same. |
---|
743 | #+nil |
---|
744 | (defun mul-qd (a b) |
---|
745 | (declare (type %quad-double a b) |
---|
746 | (optimize (speed 3))) |
---|
747 | (multiple-value-bind (a0 a1 a2 a3) |
---|
748 | (qd-parts a) |
---|
749 | (multiple-value-bind (b0 b1 b2 b3) |
---|
750 | (qd-parts b) |
---|
751 | (multiple-value-bind (p0 q0) |
---|
752 | (two-prod a0 b0) |
---|
753 | (multiple-value-bind (p1 q1) |
---|
754 | (two-prod a0 b1) |
---|
755 | (multiple-value-bind (p2 q2) |
---|
756 | (two-prod a1 b0) |
---|
757 | (multiple-value-bind (p3 q3) |
---|
758 | (two-prod a0 b2) |
---|
759 | (multiple-value-bind (p4 q4) |
---|
760 | (two-prod a1 b1) |
---|
761 | (declare (double-float q4)) |
---|
762 | #+nil |
---|
763 | (progn |
---|
764 | (format t"p0, q0 = ~a ~a~%" p0 q0) |
---|
765 | (format t"p1, q1 = ~a ~a~%" p1 q1) |
---|
766 | (format t"p2, q2 = ~a ~a~%" p2 q2) |
---|
767 | (format t"p3, q3 = ~a ~a~%" p3 q3) |
---|
768 | (format t"p4, q4 = ~a ~a~%" p4 q4)) |
---|
769 | (multiple-value-bind (p5 q5) |
---|
770 | (two-prod a2 b0) |
---|
771 | #+nil |
---|
772 | (format t"p5, q5 = ~a ~a~%" p5 q5) |
---|
773 | ;; Start accumulation |
---|
774 | (multiple-value-setq (p1 p2 q0) |
---|
775 | (three-sum p1 p2 q0)) |
---|
776 | #+nil |
---|
777 | (format t "p1 p2 q0 = ~a ~a ~a~%" p1 p2 q0) |
---|
778 | ;; six-three-sum of p2, q1, q2, p3, p4, p5 |
---|
779 | (multiple-value-setq (p2 q1 q2) |
---|
780 | (three-sum p2 q1 q2)) |
---|
781 | (multiple-value-setq (p3 p4 p5) |
---|
782 | (three-sum p3 p4 p5)) |
---|
783 | ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5) |
---|
784 | (multiple-value-bind (s0 t0) |
---|
785 | (two-sum p2 p3) |
---|
786 | (multiple-value-bind (s1 t1) |
---|
787 | (two-sum q1 p4) |
---|
788 | (declare (double-float t1)) |
---|
789 | (let ((s2 (cl:+ q2 p5))) |
---|
790 | (declare (double-float s2)) |
---|
791 | (multiple-value-bind (s1 t0) |
---|
792 | (two-sum s1 t0) |
---|
793 | (declare (double-float s1)) |
---|
794 | (incf s2 (cl:+ t0 t1)) |
---|
795 | (multiple-value-bind (p6 q6) |
---|
796 | (two-prod a0 b3) |
---|
797 | (multiple-value-bind (p7 q7) |
---|
798 | (two-prod a1 b2) |
---|
799 | (multiple-value-bind (p8 q8) |
---|
800 | (two-prod a2 b1) |
---|
801 | (multiple-value-bind (p9 q9) |
---|
802 | (two-prod a3 b0) |
---|
803 | (multiple-value-setq (q0 q3) |
---|
804 | (two-sum q0 q3)) |
---|
805 | (multiple-value-setq (q4 q5) |
---|
806 | (two-sum q4 q5)) |
---|
807 | (multiple-value-setq (p6 p7) |
---|
808 | (two-sum p6 p7)) |
---|
809 | (multiple-value-setq (p8 p9) |
---|
810 | (two-sum p8 p9)) |
---|
811 | ;; (t0,t1) = (q0,q3)+(q4,q5) |
---|
812 | (multiple-value-setq (t0 t1) |
---|
813 | (two-sum q0 q4)) |
---|
814 | (setf t1 (cl:+ q3 q5)) |
---|
815 | ;; (r0,r1) = (p6,p7)+(p8,p9) |
---|
816 | (multiple-value-bind (r0 r1) |
---|
817 | (two-sum p6 p8) |
---|
818 | (declare (double-float r1)) |
---|
819 | (incf r1 (cl:+ p7 p9)) |
---|
820 | ;; (q3,q4) = (t0,t1)+(r0,r1) |
---|
821 | (multiple-value-setq (q3 q4) |
---|
822 | (two-sum t0 r0)) |
---|
823 | (incf q4 (cl:+ t1 r1)) |
---|
824 | ;; (t0,t1)=(q3,q4)+s1 |
---|
825 | (multiple-value-setq (t0 t1) |
---|
826 | (two-sum q3 s1)) |
---|
827 | (incf t1 q4) |
---|
828 | ;; O(eps^4) terms |
---|
829 | (incf t1 |
---|
830 | (cl:+ (cl:* a1 b3) |
---|
831 | (cl:* a2 b2) |
---|
832 | (cl:* a3 b1) |
---|
833 | q6 q7 q8 q9 |
---|
834 | s2)) |
---|
835 | #+nil |
---|
836 | (format t "p0,p1,s0,t0,t1 = ~a ~a ~a ~a ~a~%" |
---|
837 | p0 p1 s0 t0 t1) |
---|
838 | (multiple-value-call #'%make-qd-d |
---|
839 | (renorm-5 p0 p1 s0 t0 t1)))))))))))))))))))) |
---|
840 | |
---|
841 | (defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
842 | "Return the square of the %QUAD-DOUBLE number A. If TARGET is also given, |
---|
843 | it is destructively modified with the result." |
---|
844 | (sqr-qd-t a target)) |
---|
845 | |
---|
846 | (defun sqr-qd-t (a target) |
---|
847 | "Square A" |
---|
848 | (declare (type %quad-double a #+oct-array target) |
---|
849 | (optimize (speed 3) |
---|
850 | (space 0)) |
---|
851 | #+(and cmu (not oct-array)) |
---|
852 | (ignore target)) |
---|
853 | (multiple-value-bind (p0 q0) |
---|
854 | (two-sqr (qd-0 a)) |
---|
855 | (multiple-value-bind (p1 q1) |
---|
856 | (two-prod (cl:* 2 (qd-0 a)) (qd-1 a)) |
---|
857 | (multiple-value-bind (p2 q2) |
---|
858 | (two-prod (cl:* 2 (qd-0 a)) (qd-2 a)) |
---|
859 | (multiple-value-bind (p3 q3) |
---|
860 | (two-sqr (qd-1 a)) |
---|
861 | (two-sum p1 q0 q0 p1) |
---|
862 | (two-sum q0 q1 q0 q1) |
---|
863 | (two-sum p2 p3 p2 p3) |
---|
864 | |
---|
865 | (let* ((s0 0d0) |
---|
866 | (t0 s0) |
---|
867 | (s1 s0) |
---|
868 | (t1 s0)) |
---|
869 | (declare (double-float s0 s1 t0 t1)) |
---|
870 | (two-sum s0 t0 q0 p2) |
---|
871 | (two-sum s1 t1 q1 p3) |
---|
872 | (two-sum s1 t0 s1 t0) |
---|
873 | (incf t0 t1) |
---|
874 | (quick-two-sum s1 t0 s1 t0) |
---|
875 | (quick-two-sum p2 t1 s0 s1) |
---|
876 | (quick-two-sum p3 q0 t1 t0) |
---|
877 | |
---|
878 | (let ((p4 (cl:* 2 (qd-0 a) (qd-3 a))) |
---|
879 | (p5 (cl:* 2 (qd-1 a) (qd-2 a)))) |
---|
880 | (declare (double-float p4 p5)) |
---|
881 | (two-sum p4 p5 p4 p5) |
---|
882 | (two-sum q2 q3 q2 q3) |
---|
883 | (two-sum t0 t1 p4 q2) |
---|
884 | |
---|
885 | (incf t1 (cl:+ p5 q3)) |
---|
886 | |
---|
887 | (two-sum p3 p4 p3 t0) |
---|
888 | (incf p4 (cl:+ q0 t1)) |
---|
889 | |
---|
890 | (multiple-value-bind (a0 a1 a2 a3) |
---|
891 | (renorm-5 p0 p1 p2 p3 p4) |
---|
892 | (%store-qd-d target a0 a1 a2 a3))))))))) |
---|
893 | |
---|
894 | |
---|
895 | (defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
896 | "Return the quotient of the two %QUAD-DOUBLE numbers A and B. |
---|
897 | If TARGET is given, it destrutively modified with the result." |
---|
898 | (div-qd-t a b target)) |
---|
899 | |
---|
900 | #+nil |
---|
901 | (defun div-qd-t (a b target) |
---|
902 | (declare (type %quad-double a b) |
---|
903 | (optimize (speed 3) |
---|
904 | (space 0)) |
---|
905 | (inline float-infinity-p) |
---|
906 | #+cmu |
---|
907 | (ignore target)) |
---|
908 | (let ((a0 (qd-0 a)) |
---|
909 | (b0 (qd-0 b))) |
---|
910 | (let* ((q0 (cl:/ a0 b0)) |
---|
911 | (r (sub-qd a (mul-qd-d b q0))) |
---|
912 | (q1 (cl:/ (qd-0 r) b0))) |
---|
913 | |
---|
914 | (when (float-infinity-p q0) |
---|
915 | (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0))) |
---|
916 | (setf r (sub-qd r (mul-qd-d b q1))) |
---|
917 | (let ((q2 (cl:/ (qd-0 r) b0))) |
---|
918 | (setf r (sub-qd r (mul-qd-d b q2))) |
---|
919 | (let ((q3 (cl:/ (qd-0 r) b0))) |
---|
920 | (multiple-value-bind (q0 q1 q2 q3) |
---|
921 | (renorm-4 q0 q1 q2 q3) |
---|
922 | (%store-qd-d target q0 q1 q2 q3))))))) |
---|
923 | |
---|
924 | (defun div-qd-t (a b target) |
---|
925 | (declare (type %quad-double a b #+oct-array target) |
---|
926 | (optimize (speed 3) |
---|
927 | (space 0)) |
---|
928 | (inline float-infinity-p) |
---|
929 | #+(and cmu (not oct-array)) |
---|
930 | (ignore target)) |
---|
931 | (let ((a0 (qd-0 a)) |
---|
932 | (b0 (qd-0 b))) |
---|
933 | (let* ((q0 (cl:/ a0 b0)) |
---|
934 | (r (%make-qd-d 0d0 0d0 0d0 0d0))) |
---|
935 | (mul-qd-d b q0 r) |
---|
936 | (sub-qd a r r) |
---|
937 | (let* ((q1 (cl:/ (qd-0 r) b0)) |
---|
938 | (temp (mul-qd-d b q1))) |
---|
939 | (when (float-infinity-p q0) |
---|
940 | (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0))) |
---|
941 | |
---|
942 | (sub-qd r temp r) |
---|
943 | (let ((q2 (cl:/ (qd-0 r) b0))) |
---|
944 | (mul-qd-d b q2 temp) |
---|
945 | (sub-qd r temp r) |
---|
946 | (let ((q3 (cl:/ (qd-0 r) b0))) |
---|
947 | (multiple-value-bind (q0 q1 q2 q3) |
---|
948 | (renorm-4 q0 q1 q2 q3) |
---|
949 | (%store-qd-d target q0 q1 q2 q3)))))))) |
---|
950 | |
---|
951 | (declaim (inline invert-qd)) |
---|
952 | |
---|
953 | (defun invert-qd (v) |
---|
954 | ;; a quartic newton iteration for 1/v |
---|
955 | ;; to invert v, start with a good guess, x. |
---|
956 | ;; let h= 1-v*x ;; h is small |
---|
957 | ;; return x+ x*(h+h^2+h^3) . compute h3 in double-float |
---|
958 | ;; enough accuracy. |
---|
959 | |
---|
960 | (let* |
---|
961 | ((x (%make-qd-d (cl:/ (qd-0 v)) 0d0 0d0 0d0)) |
---|
962 | (h (add-qd-d (neg-qd (mul-qd v x)) |
---|
963 | 1.0d0)) |
---|
964 | (h2 (mul-qd h h)) ;also use h2 for target |
---|
965 | (h3 (* (qd-0 h) (qd-0 h2)))) |
---|
966 | (declare (type %quad-double v h h2) |
---|
967 | (double-float h3)) |
---|
968 | (add-qd x |
---|
969 | (mul-qd x |
---|
970 | (add-qd h (add-qd-d h2 h3)))))) |
---|
971 | |
---|
972 | (defun div-qd-i (a b) |
---|
973 | (declare (type %quad-double a b) |
---|
974 | (optimize (speed 3) |
---|
975 | (space 0))) |
---|
976 | (mul-qd a (invert-qd b))) |
---|
977 | |
---|
978 | ;; Non-sloppy divide |
---|
979 | #+(or) |
---|
980 | (defun div-qd (a b) |
---|
981 | (declare (type %quad-double a b)) |
---|
982 | (let ((a0 (qd-0 a)) |
---|
983 | (b0 (qd-0 b))) |
---|
984 | (let* ((q0 (cl:/ a0 b0)) |
---|
985 | (r (sub-qd a (mul-qd-d b q0))) |
---|
986 | (q1 (cl:/ (qd-0 r) b0))) |
---|
987 | (setf r (sub-qd r (mul-qd-d b q1))) |
---|
988 | (let ((q2 (cl:/ (qd-0 r) b0))) |
---|
989 | (setf r (sub-qd r (mul-qd-d b q2))) |
---|
990 | (let ((q3 (cl:/ (qd-0 r) b0))) |
---|
991 | (setf r (sub-qd r (mul-qd-d b q3))) |
---|
992 | (let ((q4 (cl:/ (qd-0 r) b0))) |
---|
993 | (multiple-value-bind (q0 q1 q2 q3) |
---|
994 | (renorm-5 q0 q1 q2 q3 q4) |
---|
995 | (%make-qd-d q0 q1 q2 q3)))))))) |
---|
996 | |
---|
997 | ;; quad-double / double |
---|
998 | (defun div-qd-d (a b) |
---|
999 | "Divides the %QUAD-DOUBLE number A by the DOUBLE-FLOAT number B." |
---|
1000 | (declare (type %quad-double a) |
---|
1001 | (double-float b) |
---|
1002 | (optimize (speed 3) |
---|
1003 | (space 0)) |
---|
1004 | (inline float-infinity-p)) |
---|
1005 | ;; Compute approximate quotient using high order doubles, then |
---|
1006 | ;; correct it 3 times using the remainder. Analogous to long |
---|
1007 | ;; division. |
---|
1008 | (let ((q0 (cl:/ (qd-0 a) b))) |
---|
1009 | |
---|
1010 | (when (float-infinity-p q0) |
---|
1011 | (return-from div-qd-d (%make-qd-d q0 0d0 0d0 0d0))) |
---|
1012 | |
---|
1013 | ;; Compute remainder a - q0 * b |
---|
1014 | (multiple-value-bind (t0 t1) |
---|
1015 | (two-prod q0 b) |
---|
1016 | (let ((r #+cmu (sub-qd-dd a (kernel:make-double-double-float t0 t1)) |
---|
1017 | #-cmu (sub-qd a (make-qd-d t0 t1 0d0 0d0)))) |
---|
1018 | ;; First correction |
---|
1019 | (let ((q1 (cl:/ (qd-0 r) b))) |
---|
1020 | (multiple-value-bind (t0 t1) |
---|
1021 | (two-prod q1 b) |
---|
1022 | (setf r #+cmu (sub-qd-dd r (kernel:make-double-double-float t0 t1)) |
---|
1023 | #-cmu (sub-qd r (make-qd-d t0 t1 0d0 0d0))) |
---|
1024 | ;; Second correction |
---|
1025 | (let ((q2 (cl:/ (qd-0 r) b))) |
---|
1026 | (multiple-value-bind (t0 t1) |
---|
1027 | (two-prod q2 b) |
---|
1028 | (setf r #+cmu (sub-qd-dd r (kernel:make-double-double-float t0 t1)) |
---|
1029 | #-cmu (sub-qd r (make-qd-d t0 t1 0d0 0d0))) |
---|
1030 | ;; Final correction |
---|
1031 | (let ((q3 (cl:/ (qd-0 r) b))) |
---|
1032 | (make-qd-d q0 q1 q2 q3)))))))))) |
---|
1033 | |
---|
1034 | ;; Sloppy version |
---|
1035 | #+cmu |
---|
1036 | (defun div-qd-dd (a b) |
---|
1037 | (declare (type %quad-double a) |
---|
1038 | (double-double-float b) |
---|
1039 | (optimize (speed 3) |
---|
1040 | (space 0)) |
---|
1041 | (inline float-infinity-p)) |
---|
1042 | (let* ((q0 (cl:/ (qd-0 a) (kernel:double-double-hi b))) |
---|
1043 | (r (sub-qd-dd a (cl:* b q0)))) |
---|
1044 | (when (float-infinity-p q0) |
---|
1045 | (return-from div-qd-dd (%make-qd-d q0 0d0 0d0 0d0))) |
---|
1046 | (let ((q1 (cl:/ (qd-0 r) (kernel:double-double-hi b)))) |
---|
1047 | (setf r (sub-qd-dd r (cl:* b q1))) |
---|
1048 | (let ((q2 (cl:/ (qd-0 r) (kernel:double-double-hi b)))) |
---|
1049 | (setf r (sub-qd-dd r (cl:* b q2))) |
---|
1050 | (let ((q3 (cl:/ (qd-0 r) (kernel:double-double-hi b)))) |
---|
1051 | (make-qd-d q0 q1 q2 q3)))))) |
---|
1052 | |
---|
1053 | #+cmu |
---|
1054 | (defun make-qd-dd (a0 a1) |
---|
1055 | "Create a %QUAD-DOUBLE from two double-double-floats. This is for CMUCL, |
---|
1056 | which supports double-double floats." |
---|
1057 | (declare (double-double-float a0 a1) |
---|
1058 | (optimize (speed 3) (space 0))) |
---|
1059 | (make-qd-d (kernel:double-double-hi a0) |
---|
1060 | (kernel:double-double-lo a0) |
---|
1061 | (kernel:double-double-hi a1) |
---|
1062 | (kernel:double-double-lo a1))) |
---|
1063 | |
---|
1064 | |
---|
1065 | #-(or qd-inline (not cmu)) |
---|
1066 | (declaim (ext:end-block)) |
---|
1067 | |
---|
1068 | (defun abs-qd (a) |
---|
1069 | "Returns the absolute value of the %QUAD-DOUBLE A" |
---|
1070 | (declare (type %quad-double a)) |
---|
1071 | (if (minusp (float-sign (qd-0 a))) |
---|
1072 | (neg-qd a) |
---|
1073 | a)) |
---|
1074 | |
---|
1075 | ;; a^n for an integer n |
---|
1076 | (defun npow (a n) |
---|
1077 | "Return N'th power of A, where A is a %QUAD-DOUBLE number and N |
---|
1078 | is a fixnum." |
---|
1079 | (declare (type %quad-double a) |
---|
1080 | (fixnum n) |
---|
1081 | (optimize (speed 3) |
---|
1082 | (space 0))) |
---|
1083 | (when (= n 0) |
---|
1084 | (return-from npow (make-qd-d 1d0))) |
---|
1085 | |
---|
1086 | (let ((r a) |
---|
1087 | (s (make-qd-d 1d0)) |
---|
1088 | (abs-n (abs n))) |
---|
1089 | (declare (type (and fixnum unsigned-byte) abs-n) |
---|
1090 | (type %quad-double r s)) |
---|
1091 | (cond ((> abs-n 1) |
---|
1092 | ;; Binary exponentiation |
---|
1093 | (loop while (plusp abs-n) |
---|
1094 | do |
---|
1095 | (when (= 1 (logand abs-n 1)) |
---|
1096 | (setf s (mul-qd s r))) |
---|
1097 | (setf abs-n (ash abs-n -1)) |
---|
1098 | (when (plusp abs-n) |
---|
1099 | (setf r (sqr-qd r))))) |
---|
1100 | (t |
---|
1101 | (setf s r))) |
---|
1102 | (if (minusp n) |
---|
1103 | (div-qd (make-qd-d 1d0) s) |
---|
1104 | s))) |
---|
1105 | |
---|
1106 | ;; The n'th root of a. n is an positive integer and a should be |
---|
1107 | ;; positive too. |
---|
1108 | (defun nroot-qd (a n) |
---|
1109 | (declare (type %quad-double a) |
---|
1110 | (type (and fixnum unsigned-byte) n) |
---|
1111 | (optimize (speed 3) |
---|
1112 | (space 0))) |
---|
1113 | ;; Use Newton's iteration to solve |
---|
1114 | ;; |
---|
1115 | ;; 1/(x^n) - a = 0 |
---|
1116 | ;; |
---|
1117 | ;; The iteration becomes |
---|
1118 | ;; |
---|
1119 | ;; x' = x + x * (1 - a * x^n)/n |
---|
1120 | ;; |
---|
1121 | ;; Since Newton's iteration converges quadratically, we only need to |
---|
1122 | ;; perform it twice. |
---|
1123 | (let ((r (make-qd-d (expt (the (double-float (0d0)) (qd-0 a)) |
---|
1124 | (cl:- (cl:/ (float n 1d0))))))) |
---|
1125 | (declare (type %quad-double r)) |
---|
1126 | (flet ((term () |
---|
1127 | (div-qd-d (mul-qd r |
---|
1128 | (add-qd-d (neg-qd (mul-qd a (npow r n))) |
---|
1129 | 1d0)) |
---|
1130 | (float n 1d0)))) |
---|
1131 | (dotimes (k 3) |
---|
1132 | (setf r (add-qd r (term)))) |
---|
1133 | (div-qd (make-qd-d 1d0) r)))) |
---|
1134 | |
---|
1135 | (defun qd-< (a b) |
---|
1136 | "Returns T if A < B, where A and B are %QUAD-DOUBLE numbers." |
---|
1137 | (or (< (qd-0 a) (qd-0 b)) |
---|
1138 | (and (= (qd-0 a) (qd-0 b)) |
---|
1139 | (or (< (qd-1 a) (qd-1 b)) |
---|
1140 | (and (= (qd-1 a) (qd-1 b)) |
---|
1141 | (or (< (qd-2 a) (qd-2 b)) |
---|
1142 | (and (= (qd-2 a) (qd-2 b)) |
---|
1143 | (< (qd-3 a) (qd-3 b))))))))) |
---|
1144 | |
---|
1145 | (defun qd-> (a b) |
---|
1146 | "Returns T if A > B, where A and B are %QUAD-DOUBLE numbers." |
---|
1147 | (or (> (qd-0 a) (qd-0 b)) |
---|
1148 | (and (= (qd-0 a) (qd-0 b)) |
---|
1149 | (or (> (qd-1 a) (qd-1 b)) |
---|
1150 | (and (= (qd-1 a) (qd-1 b)) |
---|
1151 | (or (> (qd-2 a) (qd-2 b)) |
---|
1152 | (and (= (qd-2 a) (qd-2 b)) |
---|
1153 | (> (qd-3 a) (qd-3 b))))))))) |
---|
1154 | |
---|
1155 | (defun qd-<= (a b) |
---|
1156 | "Returns T if A <= B, where A and B are %QUAD-DOUBLE numbers." |
---|
1157 | (not (qd-> a b))) |
---|
1158 | |
---|
1159 | (defun qd->= (a b) |
---|
1160 | "Returns T if A >= B, where A and B are %QUAD-DOUBLE numbers." |
---|
1161 | (not (qd-< a b))) |
---|
1162 | |
---|
1163 | (defun zerop-qd (a) |
---|
1164 | "Returns T if the %QUAD-DOUBLE number A is numerically equal to 0." |
---|
1165 | (declare (type %quad-double a)) |
---|
1166 | (zerop (qd-0 a))) |
---|
1167 | |
---|
1168 | (defun onep-qd (a) |
---|
1169 | "Returns T if the %QUAD-DOUBLE number A is numerically equal to 1." |
---|
1170 | (declare (type %quad-double a)) |
---|
1171 | (and (= (qd-0 a) 1d0) |
---|
1172 | (zerop (qd-1 a)) |
---|
1173 | (zerop (qd-2 a)) |
---|
1174 | (zerop (qd-3 a)))) |
---|
1175 | |
---|
1176 | (defun plusp-qd (a) |
---|
1177 | "Returns T if the %QUAD-DOUBLE number A is strictly positive." |
---|
1178 | (declare (type %quad-double a)) |
---|
1179 | (plusp (qd-0 a))) |
---|
1180 | |
---|
1181 | (defun minusp-qd (a) |
---|
1182 | "Returns T if the %QUAD-DOUBLE number A is strictly negative." |
---|
1183 | (declare (type %quad-double a)) |
---|
1184 | (minusp (qd-0 a))) |
---|
1185 | |
---|
1186 | (defun qd-= (a b) |
---|
1187 | "Returns T if the %QUAD-DOUBLE numbers A and B are numerically equal." |
---|
1188 | (and (= (qd-0 a) (qd-0 b)) |
---|
1189 | (= (qd-1 a) (qd-1 b)) |
---|
1190 | (= (qd-2 a) (qd-2 b)) |
---|
1191 | (= (qd-3 a) (qd-3 b)))) |
---|
1192 | |
---|
1193 | |
---|
1194 | #+nil |
---|
1195 | (defun integer-decode-qd (q) |
---|
1196 | (declare (type %quad-double q)) |
---|
1197 | (multiple-value-bind (hi-int hi-exp sign) |
---|
1198 | (integer-decode-float (realpart q)) |
---|
1199 | (if (zerop (imagpart q)) |
---|
1200 | (values (ash hi-int 106) (cl:- hi-exp 106) sign) |
---|
1201 | (multiple-value-bind (lo-int lo-exp lo-sign) |
---|
1202 | (integer-decode-float (imagpart q)) |
---|
1203 | (values (cl:+ (cl:* (cl:* sign lo-sign) lo-int) |
---|
1204 | (ash hi-int (cl:- hi-exp lo-exp))) |
---|
1205 | lo-exp |
---|
1206 | sign))))) |
---|
1207 | |
---|
1208 | #+(or) |
---|
1209 | (defun integer-decode-qd (q) |
---|
1210 | (declare (type %quad-double q)) |
---|
1211 | ;; Integer decode each component and then create the appropriate |
---|
1212 | ;; integer by shifting and add all the parts together. |
---|
1213 | (multiple-value-bind (q0-int q0-exp q0-sign) |
---|
1214 | (integer-decode-float (qd-0 q)) |
---|
1215 | (multiple-value-bind (q1-int q1-exp q1-sign) |
---|
1216 | (integer-decode-float (qd-1 q)) |
---|
1217 | ;; Note: Some systems return an exponent of 0 if the number is |
---|
1218 | ;; zero. If so, everything is easier if we pretend the exponent |
---|
1219 | ;; is -1075. |
---|
1220 | (when (zerop (qd-1 q)) |
---|
1221 | (setf q1-exp -1075)) |
---|
1222 | (multiple-value-bind (q2-int q2-exp q2-sign) |
---|
1223 | (integer-decode-float (qd-2 q)) |
---|
1224 | (when (zerop (qd-2 q)) |
---|
1225 | (setf q2-exp -1075)) |
---|
1226 | (multiple-value-bind (q3-int q3-exp q3-sign) |
---|
1227 | (integer-decode-float (qd-3 q)) |
---|
1228 | (when (zerop (qd-3 q)) |
---|
1229 | (setf q3-exp -1075)) |
---|
1230 | ;; Combine all the parts together. |
---|
1231 | (values (+ (* q0-sign q3-sign q3-int) |
---|
1232 | (ash (* q0-sign q2-sign q2-int) (- q2-exp q3-exp)) |
---|
1233 | (ash (* q0-sign q1-sign q1-int) (- q1-exp q3-exp)) |
---|
1234 | (ash q0-int (- q0-exp q3-exp))) |
---|
1235 | q3-exp |
---|
1236 | q0-sign)))))) |
---|
1237 | |
---|
1238 | #+(or) |
---|
1239 | (defun integer-decode-qd (q) |
---|
1240 | (declare (type %quad-double q)) |
---|
1241 | ;; Integer decode each component and then create the appropriate |
---|
1242 | ;; integer by shifting and adding all the parts together. If any |
---|
1243 | ;; component is zero, we stop. |
---|
1244 | (multiple-value-bind (q0-int q0-exp q0-sign) |
---|
1245 | (integer-decode-float (qd-0 q)) |
---|
1246 | (if (zerop (qd-1 q)) |
---|
1247 | (values q0-int q0-exp q0-sign) |
---|
1248 | (multiple-value-bind (q1-int q1-exp q1-sign) |
---|
1249 | (integer-decode-float (qd-1 q)) |
---|
1250 | (setf q0-int (+ (ash q0-int (- q0-exp q1-exp)) |
---|
1251 | (* q1-sign q1-int))) |
---|
1252 | (if (zerop (qd-2 q)) |
---|
1253 | (values q0-int q1-exp q0-sign) |
---|
1254 | (multiple-value-bind (q2-int q2-exp q2-sign) |
---|
1255 | (integer-decode-float (qd-2 q)) |
---|
1256 | (setf q0-int (+ (ash q0-int (- q1-exp q2-exp)) |
---|
1257 | (* q2-sign q2-int))) |
---|
1258 | (if (zerop (qd-3 q)) |
---|
1259 | (values q0-int q2-exp q0-sign) |
---|
1260 | (multiple-value-bind (q3-int q3-exp q3-sign) |
---|
1261 | (integer-decode-float (qd-3 q)) |
---|
1262 | (values (+ (ash q0-int (- q2-exp q3-exp)) |
---|
1263 | (* q3-sign q3-int)) |
---|
1264 | q3-exp |
---|
1265 | q0-sign))))))))) |
---|
1266 | |
---|
1267 | (defun integer-decode-qd (q) |
---|
1268 | "Like INTEGER-DECODE-FLOAT, but for %QUAD-DOUBLE numbers. |
---|
1269 | Returns three values: |
---|
1270 | 1) an integer representation of the significand. |
---|
1271 | 2) the exponent for the power of 2 that the significand must be multiplied |
---|
1272 | by to get the actual value. |
---|
1273 | 3) -1 or 1 (i.e. the sign of the argument.)" |
---|
1274 | |
---|
1275 | (declare (type %quad-double q)) |
---|
1276 | ;; Integer decode each component and then create the appropriate |
---|
1277 | ;; integer by shifting and adding all the parts together. If any |
---|
1278 | ;; component is zero, we stop. |
---|
1279 | (multiple-value-bind (q0-int q0-exp q0-sign) |
---|
1280 | (integer-decode-float (qd-0 q)) |
---|
1281 | (when (zerop (qd-1 q)) |
---|
1282 | (return-from integer-decode-qd (values q0-int q0-exp q0-sign))) |
---|
1283 | (multiple-value-bind (q1-int q1-exp q1-sign) |
---|
1284 | (integer-decode-float (qd-1 q)) |
---|
1285 | (setf q0-int (+ (ash q0-int (- q0-exp q1-exp)) |
---|
1286 | (* q0-sign q1-sign q1-int))) |
---|
1287 | (when (zerop (qd-2 q)) |
---|
1288 | (return-from integer-decode-qd (values q0-int q1-exp q0-sign))) |
---|
1289 | (multiple-value-bind (q2-int q2-exp q2-sign) |
---|
1290 | (integer-decode-float (qd-2 q)) |
---|
1291 | (setf q0-int (+ (ash q0-int (- q1-exp q2-exp)) |
---|
1292 | (* q0-sign q2-sign q2-int))) |
---|
1293 | (when (zerop (qd-3 q)) |
---|
1294 | (return-from integer-decode-qd (values q0-int q2-exp q0-sign))) |
---|
1295 | (multiple-value-bind (q3-int q3-exp q3-sign) |
---|
1296 | (integer-decode-float (qd-3 q)) |
---|
1297 | (values (+ (ash q0-int (- q2-exp q3-exp)) |
---|
1298 | (* q0-sign q3-sign q3-int)) |
---|
1299 | q3-exp |
---|
1300 | q0-sign)))))) |
---|
1301 | |
---|
1302 | (declaim (inline scale-float-qd)) |
---|
1303 | #+(or) |
---|
1304 | (defun scale-float-qd (qd k) |
---|
1305 | (declare (type %quad-double qd) |
---|
1306 | (type fixnum k) |
---|
1307 | (optimize (speed 3) (space 0))) |
---|
1308 | ;; (space 0) to get scale-double-float inlined |
---|
1309 | (multiple-value-bind (a0 a1 a2 a3) |
---|
1310 | (qd-parts qd) |
---|
1311 | (make-qd-d (scale-float a0 k) |
---|
1312 | (scale-float a1 k) |
---|
1313 | (scale-float a2 k) |
---|
1314 | (scale-float a3 k)))) |
---|
1315 | |
---|
1316 | ;; The following method, which is faster doesn't work if QD is very |
---|
1317 | ;; large and k is very negative because we get zero as the answer, |
---|
1318 | ;; when it shouldn't be. |
---|
1319 | #+(or) |
---|
1320 | (defun scale-float-qd (qd k) |
---|
1321 | (declare (type %quad-double qd) |
---|
1322 | ;;(type (integer -1022 1022) k) |
---|
1323 | (optimize (speed 3) (space 0))) |
---|
1324 | ;; (space 0) to get scale-double-float inlined |
---|
1325 | (let ((scale (scale-float 1d0 k))) |
---|
1326 | (%make-qd-d (cl:* (qd-0 qd) scale) |
---|
1327 | (cl:* (qd-1 qd) scale) |
---|
1328 | (cl:* (qd-2 qd) scale) |
---|
1329 | (cl:* (qd-3 qd) scale)))) |
---|
1330 | |
---|
1331 | (defun scale-float-qd (qd k) |
---|
1332 | "Scale the %QUAD-DOUBLE number QD by 2^K. Just like SCALE-FLOAT" |
---|
1333 | (declare (type %quad-double qd) |
---|
1334 | (type (integer -2000 2000) k) |
---|
1335 | (optimize (speed 3) (space 0))) |
---|
1336 | ;; Split the exponent in half and do the scaling in two parts. |
---|
1337 | ;; Requires 2 multiplications, but should not prematurely return 0, |
---|
1338 | ;; and should be faster than the original version above. |
---|
1339 | (let* ((k1 (floor k 2)) |
---|
1340 | (k2 (- k k1)) |
---|
1341 | (s1 (scale-float 1d0 k1)) |
---|
1342 | (s2 (scale-float 1d0 k2))) |
---|
1343 | (with-qd-parts (q0 q1 q2 q3) |
---|
1344 | qd |
---|
1345 | (declare (double-float q0 q1 q2 q3)) |
---|
1346 | (%make-qd-d (cl:* (cl:* q0 s1) s2) |
---|
1347 | (cl:* (cl:* q1 s1) s2) |
---|
1348 | (cl:* (cl:* q2 s1) s2) |
---|
1349 | (cl:* (cl:* q3 s1) s2))))) |
---|
1350 | |
---|
1351 | (defun decode-float-qd (q) |
---|
1352 | "Like DECODE-FLOAT, but for %QUAD-DOUBLE numbers. Returns three values: |
---|
1353 | 1) a %QUAD-DOUBLE number representing the significand. This is always |
---|
1354 | between 0.5 (inclusive) and 1.0 (exclusive). |
---|
1355 | 2) an integer representing the exponent. |
---|
1356 | 3) -1.0 or 1.0 (i.e. the sign of the argument.)" |
---|
1357 | (declare (type %quad-double q)) |
---|
1358 | (multiple-value-bind (frac exp sign) |
---|
1359 | (decode-float (qd-0 q)) |
---|
1360 | (declare (ignore frac)) |
---|
1361 | ;; Got the exponent. Scale the quad-double appropriately. |
---|
1362 | (values (scale-float-qd q (- exp)) |
---|
1363 | exp |
---|
1364 | (make-qd-d sign)))) |
---|