root/qd.lisp @ 85f3e2b22f917ff35d080180ab1e2dda4476e881

Revision 85f3e2b22f917ff35d080180ab1e2dda4476e881, 39.8 KB (checked in by Raymond Toy <toy.raymond@…>, 3 years ago)

Move inline function NEG-QD-T before first use.

  • Property mode set to 100644
Line 
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
294normalizing the result from the four double-floats.  A0 is the most
295significant part of the %QUAD-DOUBLE, and A3 is the least.  The optional
296parameters 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.
313If 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.
345If 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.
414If 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.
489If 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.
496If 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
511the DOUBLE-FLOAT number B.
512
513If 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
521the %QUAD-DOUBLE number B.
522
523If 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
540the DOUBLE-FLOAT number B.
541
542If 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.
671If 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,
843it 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.
897If 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,
1056which 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
1078is 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))))
Note: See TracBrowser for help on using the browser.