source: qd.lisp @ 258cd8

cvs-ids
Last change on this file since 258cd8 was 258cd8, checked in by Raymond Toy <toy.raymond@…>, 4 years ago

Remove extra blank lines. Fix #1.

  • Property mode set to 100644
File size: 39.8 KB
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 (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
480  "Return the negative of the %QUAD-DOUBLE number A.
481If TARGET is given, TARGET is destructively modified to contain the result."
482  (neg-qd-t a target))
483
484(defun neg-qd-t (a target)
485  (declare (type %quad-double a #+oct-array target)
486           #+(and cmu (not oct-array)) (ignore target))
487  (with-qd-parts (a0 a1 a2 a3)
488      a
489    (declare (double-float a0 a1 a2 a3))
490    (%store-qd-d target (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3))))
491
492(defun sub-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
493  "Return the difference between the %QUAD-DOUBLE numbers A and B.
494If TARGET is given, TARGET is destructively modified to contain the result."
495  (declare (type %quad-double a b))
496  (add-qd-t a (neg-qd b) target))
497
498(defun sub-qd-t (a b target)
499  (add-qd-t a (neg-qd b) target))
500
501#+cmu
502(defun sub-qd-dd (a b)
503  (declare (type %quad-double a)
504           (type double-double-float b))
505  (add-qd-dd a (cl:- b)))
506
507(defun sub-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
508  "Return the difference between the %QUAD-DOUBLE number A and
509the DOUBLE-FLOAT number B.
510
511If TARGET is given, TARGET is destructively modified to contain the result."
512  (declare (type %quad-double a)
513           (type double-float b)
514           #+(and cmu (not oct-array)) (ignore target))
515  (add-qd-d a (cl:- b) #+oct-array target))
516
517(defun sub-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
518  "Return the difference between the DOUBLE-FLOAT number A and
519the %QUAD-DOUBLE number B.
520
521If TARGET is given, TARGET is destructively modified to contain the result."
522  (declare (type double-float a)
523           (type %quad-double b)
524           #+(and cmu (not oct-array)) (ignore target))
525  ;; a - b = a + (-b)
526  (add-d-qd a (neg-qd b) #+oct-array target))
527 
528
529;; Works
530;; (mul-qd-d (sqrt-qd (make-qd-dd 2w0 0w0)) 10d0) ->
531;; 14.1421356237309504880168872420969807856967187537694807317667973799q0
532;;
533;; Clisp says
534;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
535;;
536(defun mul-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
537  "Return the product of the %QUAD-DOUBLE number A and
538the DOUBLE-FLOAT number B.
539
540If TARGET is given, TARGET is destructively modified to contain the result."
541  (mul-qd-d-t a b target))
542
543(defun mul-qd-d-t (a b target)
544  "Multiply quad-double A with B"
545  (declare (type %quad-double a #+oct-array target)
546           (double-float b)
547           (optimize (speed 3)
548                     (space 0))
549           (inline float-infinity-p)
550           #+(and cmu (not oct-array)) (ignore target))
551  (multiple-value-bind (p0 q0)
552      (two-prod (qd-0 a) b)
553
554    (when (float-infinity-p p0)
555      (return-from mul-qd-d-t (%store-qd-d target p0 0d0 0d0 0d0)))
556    (multiple-value-bind (p1 q1)
557        (two-prod (qd-1 a) b)
558      (declare (double-float p1 q1))
559      (multiple-value-bind (p2 q2)
560          (two-prod (qd-2 a) b)
561        (declare (double-float p2 q2))
562        (let* ((p3 (cl:* (qd-3 a) b))
563               (s0 p0)
564               (s1 p0)
565               (s2 p0))
566          (declare (double-float s0 s1 s2 p3))
567          (two-sum s1 s2 q0 p1)
568          (three-sum s2 q1 p2 s2 q1 p2)
569          (three-sum2 q1 q2 q1 q2 p3)
570          (let ((s3 q1)
571                (s4 (cl:+ q2 p2)))
572            (multiple-value-bind (s0 s1 s2 s3)
573                (renorm-5 s0 s1 s2 s3 s4)
574              (if (zerop s0)
575                  (%store-qd-d target (float-sign p0 0d0) 0d0 0d0 0d0)
576                  (%store-qd-d target s0 s1 s2 s3)))))))))
577
578;; a0 * b0                        0
579;;      a0 * b1                   1
580;;      a1 * b0                   2
581;;           a1 * b1              3
582;;           a2 * b0              4
583;;                a2 * b1         5
584;;                a3 * b0         6
585;;                     a3 * b1    7
586
587;; Not working.
588;; (mul-qd-dd (sqrt-qd (make-qd-dd 2w0 0w0)) 10w0) ->
589;; 14.142135623730950488016887242097022172449805747901877456053837224q0
590;;
591;; But clisp says
592;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
593;;                                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
594;;
595;; Running a test program using qd (2.1.210) shows that we get the
596;; same wrong answer.
597#+(or)
598(defun mul-qd-dd (a b)
599  (declare (type %quad-double a)
600           (double-double-float b)
601           (optimize (speed 3)))
602  (multiple-value-bind (p0 q0)
603      (two-prod (qd-0 a) (kernel:double-double-hi b))
604    (multiple-value-bind (p1 q1)
605        (two-prod (qd-0 a) (kernel:double-double-lo b))
606      (multiple-value-bind (p2 q2)
607          (two-prod (qd-1 a) (kernel:double-double-hi b))
608        (multiple-value-bind (p3 q3)
609            (two-prod (qd-1 a) (kernel:double-double-lo b))
610          (multiple-value-bind (p4 q4)
611              (two-prod (qd-2 a) (kernel:double-double-hi b))
612            (format t "p0, q0 = ~A ~A~%" p0 q0)
613            (format t "p1, q1 = ~A ~A~%" p1 q1)
614            (format t "p2, q2 = ~A ~A~%" p2 q2)
615            (format t "p3, q3 = ~A ~A~%" p3 q3)
616            (format t "p4, q4 = ~A ~A~%" p4 q4)
617            (multiple-value-setq (p1 p2 q0)
618              (three-sum p1 p2 q0))
619            (format t "p1 = ~A~%" p1)
620            (format t "p2 = ~A~%" p2)
621            (format t "q0 = ~A~%" q0)
622            ;; five-three-sum
623            (multiple-value-setq (p2 p3 p4)
624              (three-sum p2 p3 p4))
625            (format t "p2 = ~A~%" p2)
626            (format t "p3 = ~A~%" p3)
627            (format t "p4 = ~A~%" p4)
628            (multiple-value-setq (q1 q2)
629              (two-sum q1 q2))
630            (multiple-value-bind (s0 t0)
631                (two-sum p2 q1)
632              (multiple-value-bind (s1 t1)
633                  (two-sum p3 q2)
634                (multiple-value-setq (s1 t0)
635                  (two-sum s1 t0))
636                (let ((s2 (cl:+ t0 t1 p4))
637                      (p2 s0)
638                      (p3 (cl:+ (cl:* (qd-2 a)
639                                (kernel:double-double-hi b))
640                             (cl:* (qd-3 a)
641                                (kernel:double-double-lo b))
642                             q3 q4)))
643                  (multiple-value-setq (p3 q0)
644                    (three-sum2 p3 q0 s1))
645                  (let ((p4 (cl:+ q0 s2)))
646                    (multiple-value-call #'%make-qd-d
647                      (renorm-5 p0 p1 p2 p3 p4))))))))))))
648
649;; a0 * b0                    0
650;;      a0 * b1               1
651;;      a1 * b0               2
652;;           a0 * b2          3
653;;           a1 * b1          4
654;;           a2 * b0          5
655;;                a0 * b3     6
656;;                a1 * b2     7
657;;                a2 * b1     8
658;;                a3 * b0     9
659
660;; Works
661;; (mul-qd (sqrt-qd (make-qd-dd 2w0 0w0)) (make-qd-dd 10w0 0w0)) ->
662;; 14.1421356237309504880168872420969807856967187537694807317667973799q0
663;;
664;; Clisp says
665;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
666
667(defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
668  "Returns the product of the %QUAD-DOUBLE numbers A and B.
669If TARGET is given, TARGET is destructively modified to contain the result."
670  (mul-qd-t a b target))
671
672(defun mul-qd-t (a b target)
673  (declare (type %quad-double a b #+oct-array target)
674           (optimize (speed 3)
675                     (space 0))
676           (inline float-infinity-p)
677           #+(and cmu (not oct-array))
678           (ignore target))
679  (with-qd-parts (a0 a1 a2 a3)
680      a
681    (declare (double-float a0 a1 a2 a3))
682    (with-qd-parts (b0 b1 b2 b3)
683        b
684      (declare (double-float b0 b1 b2 b3))
685      (multiple-value-bind (p0 q0)
686          (two-prod a0 b0)
687        #+cmu
688        (when (float-infinity-p p0)
689          (return-from mul-qd-t (%store-qd-d target p0 0d0 0d0 0d0)))
690        (multiple-value-bind (p1 q1)
691            (two-prod a0 b1)
692          (multiple-value-bind (p2 q2)
693              (two-prod a1 b0)
694            (multiple-value-bind (p3 q3)
695                (two-prod a0 b2)
696              (multiple-value-bind (p4 q4)
697                  (two-prod a1 b1)
698                (multiple-value-bind (p5 q5)
699                    (two-prod a2 b0)
700                  ;; Start accumulation
701                  (three-sum p1 p2 q0 p1 p2 q0)
702
703                  ;; six-three-sum of p2, q1, q2, p3, p4, p5
704                  (three-sum p2 q1 q2 p2 q1 q2)
705                  (three-sum p3 p4 p5 p3 p4 p5)
706                  ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5)
707                  (let* ((s0 0d0)
708                         (s1 s0)
709                         (t0 s0)
710                         (t1 s0))
711                    (declare (double-float s0 s1 t0 t1))
712                    (two-sum s0 t0 p2 p3)
713                    (two-sum s1 t1 q1 p4)
714                    (let ((s2 (cl:+ q2 p5)))
715                      (declare (double-float s2))
716                      (two-sum s1 t0 s1 t0)
717                      (incf s2 (cl:+ t0 t1))
718                      ;; O(eps^3) order terms.  This is the sloppy
719                      ;; multiplication version.  Should we use
720                      ;; the precise version?  It's significantly
721                      ;; more complex.
722                         
723                      (incf s1 (cl:+ (cl:* a0 b3)
724                                     (cl:* a1 b2)
725                                     (cl:* a2 b1)
726                                     (cl:* a3 b0)
727                                     q0 q3 q4 q5))
728                      #+nil
729                      (format t "p0,p1,s0,s1,s2 = ~a ~a ~a ~a ~a~%"
730                              p0 p1 s0 s1 s2)
731                      (multiple-value-bind (r0 r1 s0 s1)
732                          (renorm-5 p0 p1 s0 s1 s2)
733                        (if (zerop r0)
734                            (%store-qd-d target p0 0d0 0d0 0d0)
735                            (%store-qd-d target r0 r1 s0 s1))))))))))))))
736
737
738;; This is the non-sloppy version.  I think this works just fine, but
739;; since qd defaults to the sloppy multiplication version, we do the
740;; same.
741#+nil
742(defun mul-qd (a b)
743  (declare (type %quad-double a b)
744           (optimize (speed 3)))
745  (multiple-value-bind (a0 a1 a2 a3)
746      (qd-parts a)
747    (multiple-value-bind (b0 b1 b2 b3)
748        (qd-parts b)
749      (multiple-value-bind (p0 q0)
750          (two-prod a0 b0)
751        (multiple-value-bind (p1 q1)
752            (two-prod a0 b1)
753          (multiple-value-bind (p2 q2)
754              (two-prod a1 b0)
755            (multiple-value-bind (p3 q3)
756                (two-prod a0 b2)
757              (multiple-value-bind (p4 q4)
758                  (two-prod a1 b1)
759                (declare (double-float q4))
760                #+nil
761                (progn
762                  (format t"p0, q0 = ~a ~a~%" p0 q0)
763                  (format t"p1, q1 = ~a ~a~%" p1 q1)
764                  (format t"p2, q2 = ~a ~a~%" p2 q2)
765                  (format t"p3, q3 = ~a ~a~%" p3 q3)
766                  (format t"p4, q4 = ~a ~a~%" p4 q4))
767                (multiple-value-bind (p5 q5)
768                    (two-prod a2 b0)
769                  #+nil
770                  (format t"p5, q5 = ~a ~a~%" p5 q5)
771                  ;; Start accumulation
772                  (multiple-value-setq (p1 p2 q0)
773                    (three-sum p1 p2 q0))
774                  #+nil
775                  (format t "p1 p2 q0 = ~a ~a ~a~%" p1 p2 q0)
776                  ;; six-three-sum of p2, q1, q2, p3, p4, p5
777                  (multiple-value-setq (p2 q1 q2)
778                    (three-sum p2 q1 q2))
779                  (multiple-value-setq (p3 p4 p5)
780                    (three-sum p3 p4 p5))
781                  ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5)
782                  (multiple-value-bind (s0 t0)
783                      (two-sum p2 p3)
784                    (multiple-value-bind (s1 t1)
785                        (two-sum q1 p4)
786                      (declare (double-float t1))
787                      (let ((s2 (cl:+ q2 p5)))
788                        (declare (double-float s2))
789                        (multiple-value-bind (s1 t0)
790                            (two-sum s1 t0)
791                          (declare (double-float s1))
792                          (incf s2 (cl:+ t0 t1))
793                          (multiple-value-bind (p6 q6)
794                              (two-prod a0 b3)
795                            (multiple-value-bind (p7 q7)
796                                (two-prod a1 b2)
797                              (multiple-value-bind (p8 q8)
798                                  (two-prod a2 b1)
799                                (multiple-value-bind (p9 q9)
800                                    (two-prod a3 b0)
801                                  (multiple-value-setq (q0 q3)
802                                    (two-sum q0 q3))
803                                  (multiple-value-setq (q4 q5)
804                                    (two-sum q4 q5))
805                                  (multiple-value-setq (p6 p7)
806                                    (two-sum p6 p7))
807                                  (multiple-value-setq (p8 p9)
808                                    (two-sum p8 p9))
809                                  ;; (t0,t1) = (q0,q3)+(q4,q5)
810                                  (multiple-value-setq (t0 t1)
811                                    (two-sum q0 q4))
812                                  (setf t1 (cl:+ q3 q5))
813                                  ;; (r0,r1) = (p6,p7)+(p8,p9)
814                                  (multiple-value-bind (r0 r1)
815                                      (two-sum p6 p8)
816                                    (declare (double-float r1))
817                                    (incf r1 (cl:+ p7 p9))
818                                    ;; (q3,q4) = (t0,t1)+(r0,r1)
819                                    (multiple-value-setq (q3 q4)
820                                      (two-sum t0 r0))
821                                    (incf q4 (cl:+ t1 r1))
822                                    ;; (t0,t1)=(q3,q4)+s1
823                                    (multiple-value-setq (t0 t1)
824                                      (two-sum q3 s1))
825                                    (incf t1 q4)
826                                    ;; O(eps^4) terms
827                                    (incf t1
828                                          (cl:+ (cl:* a1 b3)
829                                             (cl:* a2 b2)
830                                             (cl:* a3 b1)
831                                             q6 q7 q8 q9
832                                             s2))
833                                    #+nil
834                                    (format t "p0,p1,s0,t0,t1 = ~a ~a ~a ~a ~a~%"
835                                            p0 p1 s0 t0 t1)
836                                    (multiple-value-call #'%make-qd-d
837                                      (renorm-5 p0 p1 s0 t0 t1))))))))))))))))))))
838
839(defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
840  "Return the square of the %QUAD-DOUBLE number A.  If TARGET is also given,
841it is destructively modified with the result."
842  (sqr-qd-t a target))
843
844(defun sqr-qd-t (a target)
845  "Square A"
846  (declare (type %quad-double a #+oct-array target)
847           (optimize (speed 3)
848                     (space 0))
849           #+(and cmu (not oct-array))
850           (ignore target))
851  (multiple-value-bind (p0 q0)
852      (two-sqr (qd-0 a))
853    (multiple-value-bind (p1 q1)
854        (two-prod (cl:* 2 (qd-0 a)) (qd-1 a))
855      (multiple-value-bind (p2 q2)
856          (two-prod (cl:* 2 (qd-0 a)) (qd-2 a))
857        (multiple-value-bind (p3 q3)
858            (two-sqr (qd-1 a))
859          (two-sum p1 q0 q0 p1)
860          (two-sum q0 q1 q0 q1)
861          (two-sum p2 p3 p2 p3)
862
863          (let* ((s0 0d0)
864                 (t0 s0)
865                 (s1 s0)
866                 (t1 s0))
867            (declare (double-float s0 s1 t0 t1))
868            (two-sum s0 t0 q0 p2)
869            (two-sum s1 t1 q1 p3)
870            (two-sum s1 t0 s1 t0)
871            (incf t0 t1)
872            (quick-two-sum s1 t0 s1 t0)
873            (quick-two-sum p2 t1 s0 s1)
874            (quick-two-sum p3 q0 t1 t0)
875
876            (let ((p4 (cl:* 2 (qd-0 a) (qd-3 a)))
877                  (p5 (cl:* 2 (qd-1 a) (qd-2 a))))
878              (declare (double-float p4 p5))
879              (two-sum p4 p5 p4 p5)
880              (two-sum q2 q3 q2 q3)
881              (two-sum t0 t1 p4 q2)
882
883              (incf t1 (cl:+ p5 q3))
884
885              (two-sum p3 p4 p3 t0)
886              (incf p4 (cl:+ q0 t1))
887
888              (multiple-value-bind (a0 a1 a2 a3)
889                  (renorm-5 p0 p1 p2 p3 p4)
890                (%store-qd-d target a0 a1 a2 a3)))))))))
891             
892
893(defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
894  "Return the quotient of the two %QUAD-DOUBLE numbers A and B.
895If TARGET is given, it destrutively modified with the result."
896  (div-qd-t a b target))
897
898#+nil
899(defun div-qd-t (a b target)
900  (declare (type %quad-double a b)
901           (optimize (speed 3)
902                     (space 0))
903           (inline float-infinity-p)
904           #+cmu
905           (ignore target))
906  (let ((a0 (qd-0 a))
907        (b0 (qd-0 b)))
908    (let* ((q0 (cl:/ a0 b0))
909           (r (sub-qd a (mul-qd-d b q0)))
910           (q1 (cl:/ (qd-0 r) b0)))
911
912      (when (float-infinity-p q0)
913        (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0)))
914      (setf r (sub-qd r (mul-qd-d b q1)))
915      (let ((q2 (cl:/ (qd-0 r) b0)))
916        (setf r (sub-qd r (mul-qd-d b q2)))
917        (let ((q3 (cl:/ (qd-0 r) b0)))
918          (multiple-value-bind (q0 q1 q2 q3)
919              (renorm-4 q0 q1 q2 q3)
920            (%store-qd-d target q0 q1 q2 q3)))))))
921
922(defun div-qd-t (a b target)
923  (declare (type %quad-double a b #+oct-array target)
924           (optimize (speed 3)
925                     (space 0))
926           (inline float-infinity-p)
927           #+(and cmu (not oct-array))
928           (ignore target))
929  (let ((a0 (qd-0 a))
930        (b0 (qd-0 b)))
931    (let* ((q0 (cl:/ a0 b0))
932           (r (%make-qd-d 0d0 0d0 0d0 0d0)))
933      (mul-qd-d b q0 r)
934      (sub-qd a r r)
935      (let* ((q1 (cl:/ (qd-0 r) b0))
936             (temp (mul-qd-d b q1)))
937        (when (float-infinity-p q0)
938          (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0)))
939       
940        (sub-qd r temp r)
941        (let ((q2 (cl:/ (qd-0 r) b0)))
942          (mul-qd-d b q2 temp)
943          (sub-qd r temp r)
944          (let ((q3 (cl:/ (qd-0 r) b0)))
945            (multiple-value-bind (q0 q1 q2 q3)
946                (renorm-4 q0 q1 q2 q3)
947              (%store-qd-d target q0 q1 q2 q3))))))))
948
949(declaim (inline invert-qd))
950
951(defun invert-qd (v)
952  ;; a quartic newton iteration for 1/v
953  ;; to invert v, start with a good guess, x.
954  ;; let h= 1-v*x  ;; h is small
955  ;; return x+ x*(h+h^2+h^3) . compute h3 in double-float
956  ;; enough accuracy.
957   
958  (let*
959      ((x (%make-qd-d (cl:/ (qd-0 v)) 0d0 0d0 0d0))
960       (h (add-qd-d (neg-qd (mul-qd v x))
961                    1.0d0))
962       (h2 (mul-qd h h)) ;also use h2 for target
963       (h3 (* (qd-0 h) (qd-0 h2))))
964    (declare (type %quad-double v h h2)
965             (double-float h3))
966    (add-qd x
967            (mul-qd x
968                    (add-qd h (add-qd-d h2 h3))))))
969
970(defun div-qd-i (a b)
971  (declare (type %quad-double a b)
972           (optimize (speed 3)
973                     (space 0)))
974  (mul-qd a (invert-qd b)))
975 
976;; Non-sloppy divide
977#+(or)
978(defun div-qd (a b)
979  (declare (type %quad-double a b))
980  (let ((a0 (qd-0 a))
981        (b0 (qd-0 b)))
982    (let* ((q0 (cl:/ a0 b0))
983           (r (sub-qd a (mul-qd-d b q0)))
984           (q1 (cl:/ (qd-0 r) b0)))
985      (setf r (sub-qd r (mul-qd-d b q1)))
986      (let ((q2 (cl:/ (qd-0 r) b0)))
987        (setf r (sub-qd r (mul-qd-d b q2)))
988        (let ((q3 (cl:/ (qd-0 r) b0)))
989          (setf r (sub-qd r (mul-qd-d b q3)))
990          (let ((q4 (cl:/ (qd-0 r) b0)))
991          (multiple-value-bind (q0 q1 q2 q3)
992              (renorm-5 q0 q1 q2 q3 q4)
993            (%make-qd-d q0 q1 q2 q3))))))))
994
995;; quad-double / double
996(defun div-qd-d (a b)
997  "Divides the %QUAD-DOUBLE number A by the DOUBLE-FLOAT number B."
998  (declare (type %quad-double a)
999           (double-float b)
1000           (optimize (speed 3)
1001                     (space 0))
1002           (inline float-infinity-p))
1003  ;; Compute approximate quotient using high order doubles, then
1004  ;; correct it 3 times using the remainder.  Analogous to long
1005  ;; division.
1006  (let ((q0 (cl:/ (qd-0 a) b)))
1007
1008    (when (float-infinity-p q0)
1009      (return-from div-qd-d (%make-qd-d q0 0d0 0d0 0d0)))
1010
1011    ;; Compute remainder a - q0 * b
1012    (multiple-value-bind (t0 t1)
1013        (two-prod q0 b)
1014      (let ((r #+cmu (sub-qd-dd a (kernel:make-double-double-float t0 t1))
1015               #-cmu (sub-qd a (make-qd-d t0 t1 0d0 0d0))))
1016        ;; First correction
1017        (let ((q1 (cl:/ (qd-0 r) b)))
1018          (multiple-value-bind (t0 t1)
1019              (two-prod q1 b)
1020            (setf r #+cmu (sub-qd-dd r (kernel:make-double-double-float t0 t1))
1021                    #-cmu (sub-qd r (make-qd-d t0 t1 0d0 0d0)))
1022            ;; Second correction
1023            (let ((q2 (cl:/ (qd-0 r) b)))
1024              (multiple-value-bind (t0 t1)
1025                  (two-prod q2 b)
1026                (setf r #+cmu (sub-qd-dd r (kernel:make-double-double-float t0 t1))
1027                        #-cmu (sub-qd r (make-qd-d t0 t1 0d0 0d0)))
1028                ;; Final correction
1029                (let ((q3 (cl:/ (qd-0 r) b)))
1030                  (make-qd-d q0 q1 q2 q3))))))))))
1031
1032;; Sloppy version
1033#+cmu
1034(defun div-qd-dd (a b)
1035  (declare (type %quad-double a)
1036           (double-double-float b)
1037           (optimize (speed 3)
1038                     (space 0))
1039           (inline float-infinity-p))
1040  (let* ((q0 (cl:/ (qd-0 a) (kernel:double-double-hi b)))
1041         (r (sub-qd-dd a (cl:* b q0))))
1042    (when (float-infinity-p q0)
1043      (return-from div-qd-dd (%make-qd-d q0 0d0 0d0 0d0)))
1044    (let ((q1 (cl:/ (qd-0 r) (kernel:double-double-hi b))))
1045      (setf r (sub-qd-dd r (cl:* b q1)))
1046      (let ((q2 (cl:/ (qd-0 r) (kernel:double-double-hi b))))
1047        (setf r (sub-qd-dd r (cl:* b q2)))
1048        (let ((q3 (cl:/ (qd-0 r) (kernel:double-double-hi b))))
1049          (make-qd-d q0 q1 q2 q3))))))
1050
1051#+cmu
1052(defun make-qd-dd (a0 a1)
1053  "Create a %QUAD-DOUBLE from two double-double-floats.  This is for CMUCL,
1054which supports double-double floats."
1055  (declare (double-double-float a0 a1)
1056           (optimize (speed 3) (space 0)))
1057  (make-qd-d (kernel:double-double-hi a0)
1058             (kernel:double-double-lo a0)
1059             (kernel:double-double-hi a1)
1060             (kernel:double-double-lo a1)))
1061
1062
1063#-(or qd-inline (not cmu))
1064(declaim (ext:end-block))
1065
1066(defun abs-qd (a)
1067  "Returns the absolute value of the %QUAD-DOUBLE A"
1068  (declare (type %quad-double a))
1069  (if (minusp (float-sign (qd-0 a)))
1070      (neg-qd a)
1071      a))
1072
1073;; a^n for an integer n
1074(defun npow (a n)
1075  "Return N'th power of A, where A is a %QUAD-DOUBLE number and N
1076is a fixnum."
1077  (declare (type %quad-double a)
1078           (fixnum n)
1079           (optimize (speed 3)
1080                     (space 0)))
1081  (when (= n 0)
1082    (return-from npow (make-qd-d 1d0)))
1083
1084  (let ((r a)
1085        (s (make-qd-d 1d0))
1086        (abs-n (abs n)))
1087    (declare (type (and fixnum unsigned-byte) abs-n)
1088             (type %quad-double r s))
1089    (cond ((> abs-n 1)
1090           ;; Binary exponentiation
1091           (loop while (plusp abs-n)
1092             do
1093             (when (= 1 (logand abs-n 1))
1094               (setf s (mul-qd s r)))
1095             (setf abs-n (ash abs-n -1))
1096             (when (plusp abs-n)
1097               (setf r (sqr-qd r)))))
1098          (t
1099           (setf s r)))
1100    (if (minusp n)
1101        (div-qd (make-qd-d 1d0) s)
1102        s)))
1103
1104;; The n'th root of a.  n is an positive integer and a should be
1105;; positive too.
1106(defun nroot-qd (a n)
1107  (declare (type %quad-double a)
1108           (type (and fixnum unsigned-byte) n)
1109           (optimize (speed 3)
1110                     (space 0)))
1111  ;; Use Newton's iteration to solve
1112  ;;
1113  ;; 1/(x^n) - a = 0
1114  ;;
1115  ;; The iteration becomes
1116  ;;
1117  ;; x' = x + x * (1 - a * x^n)/n
1118  ;;
1119  ;; Since Newton's iteration converges quadratically, we only need to
1120  ;; perform it twice.
1121  (let ((r (make-qd-d (expt (the (double-float (0d0)) (qd-0 a))
1122                            (cl:- (cl:/ (float n 1d0)))))))
1123    (declare (type %quad-double r))
1124    (flet ((term ()
1125             (div-qd-d (mul-qd r
1126                               (add-qd-d (neg-qd (mul-qd a (npow r n)))
1127                                         1d0))
1128                       (float n 1d0))))
1129      (dotimes (k 3)
1130        (setf r (add-qd r (term))))
1131      (div-qd (make-qd-d 1d0) r))))
1132
1133(defun qd-< (a b)
1134  "Returns T if A < B, where A and B are %QUAD-DOUBLE numbers."
1135  (or (< (qd-0 a) (qd-0 b))
1136      (and (= (qd-0 a) (qd-0 b))
1137           (or (< (qd-1 a) (qd-1 b))
1138               (and (= (qd-1 a) (qd-1 b))
1139                    (or (< (qd-2 a) (qd-2 b))
1140                        (and (= (qd-2 a) (qd-2 b))
1141                             (< (qd-3 a) (qd-3 b)))))))))
1142
1143(defun qd-> (a b)
1144  "Returns T if A > B, where A and B are %QUAD-DOUBLE numbers."
1145  (or (> (qd-0 a) (qd-0 b))
1146      (and (= (qd-0 a) (qd-0 b))
1147           (or (> (qd-1 a) (qd-1 b))
1148               (and (= (qd-1 a) (qd-1 b))
1149                    (or (> (qd-2 a) (qd-2 b))
1150                        (and (= (qd-2 a) (qd-2 b))
1151                             (> (qd-3 a) (qd-3 b)))))))))
1152
1153(defun qd-<= (a b)
1154  "Returns T if A <= B, where A and B are %QUAD-DOUBLE numbers."
1155  (not (qd-> a b)))
1156
1157(defun qd->= (a b)
1158  "Returns T if A >= B, where A and B are %QUAD-DOUBLE numbers."
1159  (not (qd-< a b)))
1160
1161(defun zerop-qd (a)
1162  "Returns T if the %QUAD-DOUBLE number A is numerically equal to 0."
1163  (declare (type %quad-double a))
1164  (zerop (qd-0 a)))
1165
1166(defun onep-qd (a)
1167  "Returns T if the %QUAD-DOUBLE number A is numerically equal to 1."
1168  (declare (type %quad-double a))
1169  (and (= (qd-0 a) 1d0)
1170       (zerop (qd-1 a))
1171       (zerop (qd-2 a))
1172       (zerop (qd-3 a))))
1173
1174(defun plusp-qd (a)
1175  "Returns T if the %QUAD-DOUBLE number A is strictly positive."
1176  (declare (type %quad-double a))
1177  (plusp (qd-0 a)))
1178         
1179(defun minusp-qd (a)
1180  "Returns T if the %QUAD-DOUBLE number A is strictly negative."
1181  (declare (type %quad-double a))
1182  (minusp (qd-0 a)))
1183
1184(defun qd-= (a b)
1185  "Returns T if the %QUAD-DOUBLE numbers A and B are numerically equal."
1186  (and (= (qd-0 a) (qd-0 b))
1187       (= (qd-1 a) (qd-1 b))
1188       (= (qd-2 a) (qd-2 b))
1189       (= (qd-3 a) (qd-3 b))))
1190
1191
1192#+nil
1193(defun integer-decode-qd (q)
1194  (declare (type %quad-double q))
1195  (multiple-value-bind (hi-int hi-exp sign)
1196      (integer-decode-float (realpart q))
1197    (if (zerop (imagpart q))
1198        (values (ash hi-int 106) (cl:- hi-exp 106) sign)
1199        (multiple-value-bind (lo-int lo-exp lo-sign)
1200            (integer-decode-float (imagpart q))
1201          (values (cl:+ (cl:* (cl:* sign lo-sign) lo-int)
1202                     (ash hi-int (cl:- hi-exp lo-exp)))
1203                  lo-exp
1204                  sign)))))
1205
1206#+(or)
1207(defun integer-decode-qd (q)
1208  (declare (type %quad-double q))
1209  ;; Integer decode each component and then create the appropriate
1210  ;; integer by shifting and add all the parts together.
1211  (multiple-value-bind (q0-int q0-exp q0-sign)
1212      (integer-decode-float (qd-0 q))
1213    (multiple-value-bind (q1-int q1-exp q1-sign)
1214        (integer-decode-float (qd-1 q))
1215      ;; Note: Some systems return an exponent of 0 if the number is
1216      ;; zero.  If so, everything is easier if we pretend the exponent
1217      ;; is -1075.
1218      (when (zerop (qd-1 q))
1219        (setf q1-exp -1075))
1220      (multiple-value-bind (q2-int q2-exp q2-sign)
1221          (integer-decode-float (qd-2 q))
1222        (when (zerop (qd-2 q))
1223          (setf q2-exp -1075))
1224        (multiple-value-bind (q3-int q3-exp q3-sign)
1225            (integer-decode-float (qd-3 q))
1226          (when (zerop (qd-3 q))
1227            (setf q3-exp -1075))
1228          ;; Combine all the parts together.
1229          (values (+ (* q0-sign q3-sign q3-int)
1230                     (ash (* q0-sign q2-sign q2-int) (- q2-exp q3-exp))
1231                     (ash (* q0-sign q1-sign q1-int) (- q1-exp q3-exp))
1232                     (ash q0-int (- q0-exp q3-exp)))
1233                  q3-exp
1234                  q0-sign))))))
1235
1236#+(or)
1237(defun integer-decode-qd (q)
1238  (declare (type %quad-double q))
1239  ;; Integer decode each component and then create the appropriate
1240  ;; integer by shifting and adding all the parts together.  If any
1241  ;; component is zero, we stop.
1242  (multiple-value-bind (q0-int q0-exp q0-sign)
1243      (integer-decode-float (qd-0 q))
1244    (if (zerop (qd-1 q))
1245        (values q0-int q0-exp q0-sign)
1246        (multiple-value-bind (q1-int q1-exp q1-sign)
1247            (integer-decode-float (qd-1 q))
1248          (setf q0-int (+ (ash q0-int (- q0-exp q1-exp))
1249                          (* q1-sign q1-int)))
1250          (if (zerop (qd-2 q))
1251              (values q0-int q1-exp q0-sign)
1252              (multiple-value-bind (q2-int q2-exp q2-sign)
1253                  (integer-decode-float (qd-2 q))
1254                (setf q0-int (+ (ash q0-int (- q1-exp q2-exp))
1255                                (* q2-sign q2-int)))
1256                (if (zerop (qd-3 q))
1257                    (values q0-int q2-exp q0-sign)
1258                    (multiple-value-bind (q3-int q3-exp q3-sign)
1259                        (integer-decode-float (qd-3 q))
1260                      (values (+ (ash q0-int (- q2-exp q3-exp))
1261                                 (* q3-sign q3-int))
1262                              q3-exp
1263                              q0-sign)))))))))
1264
1265(defun integer-decode-qd (q)
1266  "Like INTEGER-DECODE-FLOAT, but for %QUAD-DOUBLE numbers.
1267  Returns three values:
1268   1) an integer representation of the significand.
1269   2) the exponent for the power of 2 that the significand must be multiplied
1270      by to get the actual value.
1271   3) -1 or 1 (i.e. the sign of the argument.)"
1272
1273  (declare (type %quad-double q))
1274  ;; Integer decode each component and then create the appropriate
1275  ;; integer by shifting and adding all the parts together.  If any
1276  ;; component is zero, we stop.
1277  (multiple-value-bind (q0-int q0-exp q0-sign)
1278      (integer-decode-float (qd-0 q))
1279    (when (zerop (qd-1 q))
1280      (return-from integer-decode-qd (values q0-int q0-exp q0-sign)))
1281    (multiple-value-bind (q1-int q1-exp q1-sign)
1282        (integer-decode-float (qd-1 q))
1283      (setf q0-int (+ (ash q0-int (- q0-exp q1-exp))
1284                      (* q0-sign q1-sign q1-int)))
1285      (when (zerop (qd-2 q))
1286        (return-from integer-decode-qd (values q0-int q1-exp q0-sign)))
1287      (multiple-value-bind (q2-int q2-exp q2-sign)
1288          (integer-decode-float (qd-2 q))
1289        (setf q0-int (+ (ash q0-int (- q1-exp q2-exp))
1290                        (* q0-sign q2-sign q2-int)))
1291        (when (zerop (qd-3 q))
1292          (return-from integer-decode-qd (values q0-int q2-exp q0-sign)))
1293        (multiple-value-bind (q3-int q3-exp q3-sign)
1294            (integer-decode-float (qd-3 q))
1295          (values (+ (ash q0-int (- q2-exp q3-exp))
1296                     (* q0-sign q3-sign q3-int))
1297                  q3-exp
1298                  q0-sign))))))
1299
1300(declaim (inline scale-float-qd))
1301#+(or)
1302(defun scale-float-qd (qd k)
1303  (declare (type %quad-double qd)
1304           (type fixnum k)
1305           (optimize (speed 3) (space 0)))
1306  ;; (space 0) to get scale-double-float inlined
1307  (multiple-value-bind (a0 a1 a2 a3)
1308      (qd-parts qd)
1309    (make-qd-d (scale-float a0 k)
1310               (scale-float a1 k)
1311               (scale-float a2 k)
1312               (scale-float a3 k))))
1313
1314;; The following method, which is faster doesn't work if QD is very
1315;; large and k is very negative because we get zero as the answer,
1316;; when it shouldn't be.
1317#+(or)
1318(defun scale-float-qd (qd k)
1319  (declare (type %quad-double qd)
1320           ;;(type (integer -1022 1022) k)
1321           (optimize (speed 3) (space 0)))
1322  ;; (space 0) to get scale-double-float inlined
1323  (let ((scale (scale-float 1d0 k)))
1324    (%make-qd-d (cl:* (qd-0 qd) scale)
1325                (cl:* (qd-1 qd) scale)
1326                (cl:* (qd-2 qd) scale)
1327                (cl:* (qd-3 qd) scale))))
1328
1329(defun scale-float-qd (qd k)
1330  "Scale the %QUAD-DOUBLE number QD by 2^K.  Just like SCALE-FLOAT"
1331  (declare (type %quad-double qd)
1332           (type (integer -2000 2000) k)
1333           (optimize (speed 3) (space 0)))
1334  ;; Split the exponent in half and do the scaling in two parts.
1335  ;; Requires 2 multiplications, but should not prematurely return 0,
1336  ;; and should be faster than the original version above.
1337  (let* ((k1 (floor k 2))
1338         (k2 (- k k1))
1339         (s1 (scale-float 1d0 k1))
1340         (s2 (scale-float 1d0 k2)))
1341    (with-qd-parts (q0 q1 q2 q3)
1342        qd
1343      (declare (double-float q0 q1 q2 q3))
1344      (%make-qd-d (cl:* (cl:* q0 s1) s2)
1345                  (cl:* (cl:* q1 s1) s2)
1346                  (cl:* (cl:* q2 s1) s2)
1347                  (cl:* (cl:* q3 s1) s2)))))
1348
1349(defun decode-float-qd (q)
1350  "Like DECODE-FLOAT, but for %QUAD-DOUBLE numbers.  Returns three values:
1351   1) a %QUAD-DOUBLE number representing the significand.  This is always
1352      between 0.5 (inclusive) and 1.0 (exclusive).
1353   2) an integer representing the exponent.
1354   3) -1.0 or 1.0 (i.e. the sign of the argument.)"
1355  (declare (type %quad-double q))
1356  (multiple-value-bind (frac exp sign)
1357      (decode-float (qd-0 q))
1358    (declare (ignore frac))
1359    ;; Got the exponent.  Scale the quad-double appropriately.
1360    (values (scale-float-qd q (- exp))
1361            exp
1362            (make-qd-d sign))))
Note: See TracBrowser for help on using the repository browser.