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

;;; This file contains the core routines for basic arithmetic
;;; operations on a %quad-double.  This includes addition,
;;; subtraction, multiplication, division, negation. and absolute
;;; value.  Some additional helper functions are included such as
;;; raising to an integer power. and the n'th root of a (non-negative)
;;; %quad-double.  The basic comparison operators are included, and
;;; some simple tests for zerop, onep, plusp, and minusp. 
;;;
34
35
36
;;; The basic algorithms are based on Yozo Hida's double-double
;;; implementation.  However, some were copied from CMUCL and modified
;;; to support quad-doubles.
toy's avatar
toy committed
37

toy's avatar
toy committed
38

39
(in-package #:octi)
toy's avatar
toy committed
40

41
#+cmu
42
(eval-when (:compile-toplevel)
Raymond Toy's avatar
Raymond Toy committed
43
  (setf ext:*inline-expansion-limit* 1600))
44

45
;; All of the following functions should be inline.
46
;;(declaim (inline three-sum2))
toy's avatar
toy committed
47
48

;; Internal routines for implementing quad-double.
49
50
51
52
53
54
55

#+cmu
(defmacro two-sum (s e x y)
  `(multiple-value-setq (,s ,e)
    (c::two-sum ,x ,y)))


56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
(defmacro three-sum (s0 s1 s2 a b c)
  (let ((aa (gensym))
	(bb (gensym))
	(cc (gensym))
	(t1 (gensym))
	(t2 (gensym))
	(t3 (gensym)))
    `(let* ((,aa ,a)
	    (,bb ,b)
	    (,cc ,c)
	    (,t1 0d0)
	    (,t2 ,t1)
	    (,t3 ,t1))
      (declare (double-float ,t1 ,t2 ,t3))
      (two-sum ,t1 ,t2 ,aa ,bb)
      (two-sum ,s0 ,t3 ,cc ,t1)
      (two-sum ,s1 ,s2 ,t2 ,t3))))
      
	  
75

76
#+nil
77
78
79
80
81
82
83
84
(defun three-sum2 (a b c)
  (declare (double-float a b c)
	   (optimize (speed 3)))
  (let* ((t1 0d0)
	 (t2 t1)
	 (t3 t1))
    (two-sum t1 t2 a b)
    (two-sum a t3 c t1)
85
    (values a (cl:+ t2 t3))))
86

87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
(defmacro three-sum2 (s0 s1 a b c)
  (let ((aa (gensym))
	(bb (gensym))
	(cc (gensym))
	(t1 (gensym))
	(t2 (gensym))
	(t3 (gensym)))
    `(let* ((,aa ,a)
	    (,bb ,b)
	    (,cc ,c)
	    (,t1 0d0)
	    (,t2 ,t1)
	    (,t3 ,t1))
       (declare (double-float ,t1 ,t2 ,t3))
       (two-sum ,t1 ,t2 ,aa ,bb)
       (two-sum ,s0 ,t3 ,cc ,t1)
       (setf ,s1 (+ ,t2 ,t3)))))
104

toy's avatar
toy committed
105
106
107
108
109
110
111
112
;; Not needed????
#+nil
(declaim (inline quick-three-accum))
#+nil
(defun quick-three-accum (a b c)
  (declare (double-float a b c)
	   (optimize (speed 3) (space 0)))
  (multiple-value-bind (s b)
113
      (two-sum b c)
toy's avatar
toy committed
114
    (multiple-value-bind (s a)
115
	(two-sum a s)
toy's avatar
toy committed
116
117
118
      (let ((za (/= a 0))
	    (zb (/= b 0)))
	(when (and za zb)
119
	  (return-from quick-three-accum (values (cl:+ s 0d0) (cl:+ a 0d0) (cl:+ b 0d0))))
toy's avatar
toy committed
120
121
122
123
124
125
126
127
128
	(when (not za)
	  (setf a s)
	  (setf s 0d0))
	(when (not zb)
	  (setf b a)
	  (setf a s))
	(values 0d0 a b)))))


129
130
131
132
133
134
;; These functions are quite short, so we inline them to minimize
;; consing.
(declaim (inline make-qd-d
		 add-d-qd
		 add-dd-qd
		 neg-qd
135
		 neg-qd-t
136
137
138
139
		 sub-qd
		 sub-qd-dd
		 sub-qd-d
		 sub-d-qd
140
		 #+cmu
141
142
143
		 make-qd-dd
		 abs-qd
		 qd->
144
145
146
		 qd-<
		 qd->=
		 qd-<=
147
148
149
150
151
152
153
154
155
156
		 zerop-qd
		 onep-qd
		 plusp-qd
		 minusp-qd
		 qd-=
		 scale-float-qd))

;; Should these functions be inline?  The QD C++ source has these as
;; inline functions, but these are fairly large functions.  However,
;; inlining makes quite a big difference in speed and consing.
157
#+cmu
158
(declaim (#+qd-inline inline
Raymond Toy's avatar
Raymond Toy committed
159
	 #-qd-inline ext:maybe-inline
160
161
162
	 renorm-4
	 renorm-5
	 add-qd-d
163
	 add-qd-d-t
164
165
	 add-qd-dd
	 add-qd
166
	 add-qd-t
167
	 mul-qd-d
168
	 mul-qd-d-t
169
170
	 mul-qd-dd
	 mul-qd
171
	 mul-qd-t
172
	 sqr-qd
173
	 sqr-qd-t
174
	 div-qd
175
	 div-qd-t
176
177
178
	 div-qd-d
	 div-qd-dd))

179
180
181
182
183
#+cmu
(defmacro quick-two-sum (s e x y)
  `(multiple-value-setq (,s ,e)
    (c::quick-two-sum ,x ,y)))

184
#-(or qd-inline (not cmu))
185
(declaim (ext:start-block renorm-4 renorm-5
toy's avatar
toy committed
186
			  make-qd-d
toy's avatar
toy committed
187
188
			  add-qd-d add-d-qd add-qd-dd
			  add-dd-qd
189
			  add-qd add-qd-t
190
			  add-qd-d-t
toy's avatar
toy committed
191
			  neg-qd
192
			  sub-qd sub-qd-dd sub-qd-d sub-d-qd
193
194
195
			  mul-qd-d mul-qd-dd
			  mul-qd
			  mul-qd-t
196
			  mul-qd-d-t
toy's avatar
toy committed
197
			  sqr-qd
198
			  sqr-qd-t
toy's avatar
toy committed
199
			  div-qd div-qd-d div-qd-dd
200
			  div-qd-t
toy's avatar
toy committed
201
			  make-qd-dd
202
			  ))
toy's avatar
toy committed
203

toy's avatar
toy committed
204
#+(or)
toy's avatar
toy committed
205
206
207
208
(defun quick-renorm (c0 c1 c2 c3 c4)
  (declare (double-float c0 c1 c2 c3 c4)
	   (optimize (speed 3)))
  (multiple-value-bind (s t3)
209
      (quick-two-sum c3 c4)
toy's avatar
toy committed
210
    (multiple-value-bind (s t2)
211
	(quick-two-sum c2 s)
toy's avatar
toy committed
212
      (multiple-value-bind (s t1)
213
	  (quick-two-sum c1 s)
toy's avatar
toy committed
214
	(multiple-value-bind (c0 t0)
215
	    (quick-two-sum c0 s)
toy's avatar
toy committed
216
	  (multiple-value-bind (s t2)
217
	      (quick-two-sum t2 t3)
toy's avatar
toy committed
218
	    (multiple-value-bind (s t1)
219
		(quick-two-sum t1 s)
toy's avatar
toy committed
220
	      (multiple-value-bind (c1 t0)
221
		  (quick-two-sum t0 s)
toy's avatar
toy committed
222
		(multiple-value-bind (s t1)
223
		    (quick-two-sum t1 t2)
toy's avatar
toy committed
224
		  (multiple-value-bind (c2 t0)
225
		      (quick-two-sum t0 s)
226
		    (values c0 c1 c2 (cl:+ t0 t1))))))))))))
toy's avatar
toy committed
227
228
229

(defun renorm-4 (c0 c1 c2 c3)
  (declare (double-float c0 c1 c2 c3)
230
	   (optimize (speed 3) (safety #-allegro 0 #+allegro 1) (debug 0)))
231
232
233
  (let ((s0 0d0)
	(s1 0d0)
	(s2 0d0)
toy's avatar
toy committed
234
	(s3 0d0))
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
    (declare (double-float s0 s1 s2 s3))
    (quick-two-sum s0 c3 c2 c3)
    (quick-two-sum s0 c2 c1 s0)
    (quick-two-sum c0 c1 c0 s0)
    (setf s0 c0)
    (setf s1 c1)
    (cond ((/= s1 0)
	   (quick-two-sum s1 s2 s1 c2)
	   (if (/= s2 0)
	       (quick-two-sum s2 s3 s2 c3)
	       (quick-two-sum s1 s2 s1 c3)))
	  (t
	   (quick-two-sum s0 s1 s0 c2)
	   (if (/= s1 0)
	       (quick-two-sum s1 s2 s1 c3)
	       (quick-two-sum s0 s1 s0 c3))))
    (values s0 s1 s2 s3)))
  
toy's avatar
toy committed
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
257
258
  (let ((s0 0d0)
	(s1 0d0)
	(s2 0d0)
toy's avatar
toy committed
259
	(s3 0d0))
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
    (declare (double-float s0 s1 s2 s3))
    (quick-two-sum s0 c4 c3 c4)
    (quick-two-sum s0 c3 c2 s0)
    (quick-two-sum s0 c2 c1 s0)
    (quick-two-sum c0 c1 c0 s0)
    (quick-two-sum s0 s1 c0 c1)
    (cond ((/= s1 0)
	   (quick-two-sum s1 s2 s1 c2)
	   (cond ((/= s2 0)
		  (quick-two-sum s2 s3 s2 c3)
		  (if (/= s3 0)
		      (incf s3 c4)
		      (incf s2 c4)))
		 (t
		  (quick-two-sum s1 s2 s1 c3)
		  (if (/= s2 0)
		      (quick-two-sum s2 s3 s2 c4)
		      (quick-two-sum s1 s2 s1 c4)))))
	  (t
	   (quick-two-sum s0 s1 s0 c2)
	   (cond ((/= s1 0)
		  (quick-two-sum s1 s2 s1 c3)
		  (if (/= s2 0)
		      (quick-two-sum s2 s3 s2 c4)
		      (quick-two-sum s1 s2 s1 c4)))
		 (t
		  (quick-two-sum s0 s1 s0 c3)
		  (if (/= s1 0)
		      (quick-two-sum s1 s2 s1 c4)
		      (quick-two-sum s0 s1 s0 c4))))))
    (values s0 s1 s2 s3)))
toy's avatar
toy committed
291

toy's avatar
toy committed
292
(defun make-qd-d (a0 &optional (a1 0d0 a1-p) (a2 0d0) (a3 0d0))
Raymond Toy's avatar
Raymond Toy committed
293
294
295
296
  "Create a %QUAD-DOUBLE from four double-floats, appropriately
normalizing the result from the four double-floats.  A0 is the most
significant part of the %QUAD-DOUBLE, and A3 is the least.  The optional
parameters default to 0.
toy's avatar
toy committed
297
"
toy's avatar
toy committed
298
299
300
301
302
303
304
305
306
  (declare (double-float a0 a1 a2 a3)
	   (optimize (speed 3)
		     #+cmu
		     (ext:inhibit-warnings 3)))
  (if a1-p
      (multiple-value-bind (s0 s1 s2 s3)
	  (renorm-4 a0 a1 a2 a3)
	(%make-qd-d s0 s1 s2 s3))
      (%make-qd-d a0 0d0 0d0 0d0)))
toy's avatar
toy committed
307

toy's avatar
toy committed
308
309
310
;;;; Addition

;; Quad-double + double
311
(defun add-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
312
  "Return the sum of the %QUAD-DOUBLE A and the DOUBLE-FLOAT B.
Raymond Toy's avatar
Raymond Toy committed
313
If TARGET is given, TARGET is destructively modified to contain the result."
314
315
316
  (add-qd-d-t a b target))
  
(defun add-qd-d-t (a b target)
toy's avatar
toy committed
317
  "Add a quad-double A and a double-float B"
318
  (declare (type %quad-double a #+oct-array target)
toy's avatar
toy committed
319
	   (double-float b)
320
	   (optimize (speed 3)
Raymond Toy's avatar
Raymond Toy committed
321
		     (space 0))
322
	   (inline float-infinity-p)
323
	   #+(and cmu (not oct-array)) (ignore target))
324
325
326
327
328
329
330
  (let* ((c0 0d0)
	 (e c0)
	 (c1 c0)
	 (c2 c0)
	 (c3 c0))
    (declare (double-float e c0 c1 c2 c3))
    (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
335
336
337
338
339
    (two-sum c1 e (qd-1 a) e)
    (two-sum c2 e (qd-2 a) e)
    (two-sum c3 e (qd-3 a) e)
    (multiple-value-bind (r0 r1 r2 r3)
	(renorm-5 c0 c1 c2 c3 e)
      (if (and (zerop (qd-0 a)) (zerop b))
340
341
	  (%store-qd-d target c0 0d0 0d0 0d0)
	  (%store-qd-d target r0 r1 r2 r3)))))
toy's avatar
toy committed
342

343
(defun add-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
344
345
  "Return the sum of the DOUBLE-FLOAT A and the %QUAD-DOUBLE B.
If TARGET is given, TARGET is destructively modified to contain the result."
toy's avatar
toy committed
346
347
  (declare (double-float a)
	   (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))
toy's avatar
toy committed
351

352
#+cmu
toy's avatar
toy committed
353
354
355
356
(defun add-qd-dd (a b)
  "Add a quad-double A and a double-double B"
  (declare (type %quad-double a)
	   (double-double-float b)
357
358
	   (optimize (speed 3)
		     (space 0)))
359
360
361
362
  (let* ((s0 0d0)
	 (t0 s0)
	 (s1 s0)
	 (t1 s0)
363
	 (s2 s0)
364
	 (s3 s0))
365
    (declare (double-float s0 s1 s3 t0 t1 s2))
366
367
368
    (two-sum s0 t1 (qd-0 a) (kernel:double-double-hi b))
    (two-sum s1 t1 (qd-1 a) (kernel:double-double-lo b))
    (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
372
373
374
375
    (incf t0 t1)
    (multiple-value-bind (a0 a1 a2 a3)
	(renorm-5 s0 s1 s2 s3 t0)
      (declare (double-float a0 a1 a2 a3))
      (%make-qd-d a0 a1 a2 a3))))
toy's avatar
toy committed
376

377
#+cmu
toy's avatar
toy committed
378
379
380
(defun add-dd-qd (a b)
  (declare (double-double-float a)
	   (type %quad-double b)
381
382
	   (optimize (speed 3)
		     (space 0)))
toy's avatar
toy committed
383
384
  (add-qd-dd b a))

toy's avatar
toy committed
385

toy's avatar
toy committed
386
387
#+(or)
(defun add-qd-1 (a b)
toy's avatar
toy committed
388
389
390
  (declare (type %quad-double a b)
	   (optimize (speed 3)))
  (multiple-value-bind (s0 t0)
391
      (two-sum (qd-0 a) (qd-0 b))
toy's avatar
toy committed
392
    (multiple-value-bind (s1 t1)
393
	(two-sum (qd-1 a) (qd-1 b))
toy's avatar
toy committed
394
      (multiple-value-bind (s2 t2)
395
	  (two-sum (qd-2 a) (qd-2 b))
toy's avatar
toy committed
396
	(multiple-value-bind (s3 t3)
397
	    (two-sum (qd-3 a) (qd-3 b))
toy's avatar
toy committed
398
	  (multiple-value-bind (s1 t0)
399
	      (two-sum s1 t0)
toy's avatar
toy committed
400
401
402
403
	    (multiple-value-bind (s2 t0 t1)
		(three-sum s2 t0 t1)
	      (multiple-value-bind (s3 t0)
		  (three-sum2 s3 t0 t2)
404
		(let ((t0 (cl:+ t0 t1 t3)))
toy's avatar
toy committed
405
406
		  (multiple-value-call #'%make-qd-d
		    (renorm-5 s0 s1 s2 s3 t0)))))))))))
toy's avatar
toy committed
407
408
409
410
411

;; Same as above, except that everything is expanded out for compilers
;; which don't do a very good job with dataflow.  CMUCL is one of
;; those compilers.

412
413

(defun add-qd-t (a b target)
414
  (declare (type %quad-double a b #+oct-array target)
415
	   (optimize (speed 3)
416
		     (space 0))
417
	   #+(and cmu (not oct-array))
418
	   (ignore target))
toy's avatar
toy committed
419
420
421
422
  ;; This is the version that is NOT IEEE.  Should we use the IEEE
  ;; version?  It's quite a bit more complicated.
  ;;
  ;; In addition, this is reorganized to minimize data dependency.
Raymond Toy's avatar
Raymond Toy committed
423
424
425
426
427
428
  (with-qd-parts (a0 a1 a2 a3)
      a
    (declare (double-float a0 a1 a2 a3))
    (with-qd-parts (b0 b1 b2 b3)
	b
      (declare (double-float b0 b1 b2 b3))
429
430
431
432
      (let ((s0 (cl:+ a0 b0))
	    (s1 (cl:+ a1 b1))
	    (s2 (cl:+ a2 b2))
	    (s3 (cl:+ a3 b3)))
Raymond Toy's avatar
Raymond Toy committed
433
	(declare (double-float s0 s1 s2 s3)
434
435
		 (inline float-infinity-p))

436
	(when (float-infinity-p s0)
437
	  (return-from add-qd-t (%store-qd-d target s0 0d0 0d0 0d0)))
438
439
440
441
	(let ((v0 (cl:- s0 a0))
	      (v1 (cl:- s1 a1))
	      (v2 (cl:- s2 a2))
	      (v3 (cl:- s3 a3)))
Raymond Toy's avatar
Raymond Toy committed
442
	  (declare (double-float v0 v1 v2 v3))
443
444
445
446
	  (let ((u0 (cl:- s0 v0))
		(u1 (cl:- s1 v1))
		(u2 (cl:- s2 v2))
		(u3 (cl:- s3 v3)))
Raymond Toy's avatar
Raymond Toy committed
447
	    (declare (double-float u0 u1 u2 u3))
448
449
450
451
	    (let ((w0 (cl:- a0 u0))
		  (w1 (cl:- a1 u1))
		  (w2 (cl:- a2 u2))
		  (w3 (cl:- a3 u3)))
Raymond Toy's avatar
Raymond Toy committed
452
	      (declare (double-float w0 w1 w2 w3))
453
454
455
456
	      (let ((u0 (cl:- b0 v0))
		    (u1 (cl:- b1 v1))
		    (u2 (cl:- b2 v2))
		    (u3 (cl:- b3 v3)))
Raymond Toy's avatar
Raymond Toy committed
457
		(declare (double-float u0 u1 u2 u3))
458
459
460
461
		(let ((t0 (cl:+ w0 u0))
		      (t1 (cl:+ w1 u1))
		      (t2 (cl:+ w2 u2))
		      (t3 (cl:+ w3 u3)))
462
		  (declare (double-float t0 t1 t2 t3))
463
		  (two-sum s1 t0 s1 t0)
464
		  (three-sum s2 t0 t1 s2 t0 t1)
465
466
		  (three-sum2 s3 t0 s3 t0 t2)
		  (setf t0 (cl:+ t0 t1 t3))
Raymond Toy's avatar
Raymond Toy committed
467
468
469
470
		  ;; Renormalize
		  (multiple-value-setq (s0 s1 s2 s3)
		    (renorm-5 s0 s1 s2 s3 t0))
		  (if (and (zerop a0) (zerop b0))
471
472
473
		      (%store-qd-d target (+ a0 b0) 0d0 0d0 0d0)
		      (%store-qd-d target s0 s1 s2 s3)))))))))))

474
475
476
477
478
479
(defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
  "Return the sum of the %QUAD-DOUBLE numbers A and B.
If TARGET is given, TARGET is destructively modified to contain the result."
  (add-qd-t a b target))


480
(defun neg-qd-t (a target)
481
  (declare (type %quad-double a #+oct-array target)
482
	   #+(and cmu (not oct-array)) (ignore target))
Raymond Toy's avatar
Raymond Toy committed
483
484
485
  (with-qd-parts (a0 a1 a2 a3)
      a
    (declare (double-float a0 a1 a2 a3))
486
    (%store-qd-d target (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3))))
toy's avatar
toy committed
487

488
489
490
491
492
493
494
(defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
  "Return the negative of the %QUAD-DOUBLE number A.
If TARGET is given, TARGET is destructively modified to contain the result."
  (neg-qd-t a target))



495
(defun sub-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
496
497
  "Return the difference between the %QUAD-DOUBLE numbers A and B.
If TARGET is given, TARGET is destructively modified to contain the result."
toy's avatar
toy committed
498
  (declare (type %quad-double a b))
499
  (add-qd-t a (neg-qd b) target))
toy's avatar
toy committed
500

501
502
503
(defun sub-qd-t (a b target)
  (add-qd-t a (neg-qd b) target))

toy's avatar
toy committed
504
#+cmu
505
506
507
(defun sub-qd-dd (a b)
  (declare (type %quad-double a)
	   (type double-double-float b))
508
  (add-qd-dd a (cl:- b)))
509

510
(defun sub-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
511
512
513
514
  "Return the difference between the %QUAD-DOUBLE number A and
the DOUBLE-FLOAT number B.

If TARGET is given, TARGET is destructively modified to contain the result."
515
  (declare (type %quad-double a)
516
	   (type double-float b)
517
	   #+(and cmu (not oct-array)) (ignore target))
518
  (add-qd-d a (cl:- b) #+oct-array target))
519

520
(defun sub-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
521
522
523
524
  "Return the difference between the DOUBLE-FLOAT number A and
the %QUAD-DOUBLE number B.

If TARGET is given, TARGET is destructively modified to contain the result."
525
  (declare (type double-float a)
526
	   (type %quad-double b)
527
	   #+(and cmu (not oct-array)) (ignore target))
toy's avatar
toy committed
528
  ;; a - b = a + (-b)
529
  (add-d-qd a (neg-qd b) #+oct-array target))
530
  
toy's avatar
toy committed
531

toy's avatar
toy committed
532
533
534
535
536
537
538
;; Works
;; (mul-qd-d (sqrt-qd (make-qd-dd 2w0 0w0)) 10d0) ->
;; 14.1421356237309504880168872420969807856967187537694807317667973799q0
;;
;; Clisp says
;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
;;
539
(defun mul-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
Raymond Toy's avatar
Raymond Toy committed
540
541
542
543
  "Return the product of the %QUAD-DOUBLE number A and
the DOUBLE-FLOAT number B.

If TARGET is given, TARGET is destructively modified to contain the result."
544
545
546
  (mul-qd-d-t a b target))

(defun mul-qd-d-t (a b target)
toy's avatar
toy committed
547
  "Multiply quad-double A with B"
548
  (declare (type %quad-double a #+oct-array target)
toy's avatar
toy committed
549
	   (double-float b)
550
	   (optimize (speed 3)
Raymond Toy's avatar
Raymond Toy committed
551
		     (space 0))
552
	   (inline float-infinity-p)
553
	   #+(and cmu (not oct-array)) (ignore target))
toy's avatar
toy committed
554
  (multiple-value-bind (p0 q0)
555
      (two-prod (qd-0 a) b)
556

557
    (when (float-infinity-p p0)
558
      (return-from mul-qd-d-t (%store-qd-d target p0 0d0 0d0 0d0)))
toy's avatar
toy committed
559
    (multiple-value-bind (p1 q1)
560
	(two-prod (qd-1 a) b)
561
      (declare (double-float p1 q1))
toy's avatar
toy committed
562
      (multiple-value-bind (p2 q2)
563
	  (two-prod (qd-2 a) b)
564
	(declare (double-float p2 q2))
565
566
567
568
569
570
	(let* ((p3 (cl:* (qd-3 a) b))
	       (s0 p0)
	       (s1 p0)
	       (s2 p0))
	  (declare (double-float s0 s1 s2 p3))
	  (two-sum s1 s2 q0 p1)
571
	  (three-sum s2 q1 p2 s2 q1 p2)
572
573
574
575
576
577
	  (three-sum2 q1 q2 q1 q2 p3)
	  (let ((s3 q1)
		(s4 (cl:+ q2 p2)))
	    (multiple-value-bind (s0 s1 s2 s3)
		(renorm-5 s0 s1 s2 s3 s4)
	      (if (zerop s0)
578
579
		  (%store-qd-d target (float-sign p0 0d0) 0d0 0d0 0d0)
		  (%store-qd-d target s0 s1 s2 s3)))))))))
toy's avatar
toy committed
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599

;; a0 * b0                        0
;;      a0 * b1                   1
;;      a1 * b0                   2
;;           a1 * b1              3
;;           a2 * b0              4
;;                a2 * b1         5
;;                a3 * b0         6
;;                     a3 * b1    7

;; Not working.
;; (mul-qd-dd (sqrt-qd (make-qd-dd 2w0 0w0)) 10w0) ->
;; 14.142135623730950488016887242097022172449805747901877456053837224q0
;;
;; But clisp says
;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
;;                                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
;;
;; Running a test program using qd (2.1.210) shows that we get the
;; same wrong answer.
600
#+(or)
toy's avatar
toy committed
601
602
603
604
605
(defun mul-qd-dd (a b)
  (declare (type %quad-double a)
	   (double-double-float b)
	   (optimize (speed 3)))
  (multiple-value-bind (p0 q0)
606
      (two-prod (qd-0 a) (kernel:double-double-hi b))
toy's avatar
toy committed
607
    (multiple-value-bind (p1 q1)
608
	(two-prod (qd-0 a) (kernel:double-double-lo b))
toy's avatar
toy committed
609
      (multiple-value-bind (p2 q2)
610
	  (two-prod (qd-1 a) (kernel:double-double-hi b))
toy's avatar
toy committed
611
	(multiple-value-bind (p3 q3)
612
	    (two-prod (qd-1 a) (kernel:double-double-lo b))
toy's avatar
toy committed
613
	  (multiple-value-bind (p4 q4)
614
	      (two-prod (qd-2 a) (kernel:double-double-hi b))
toy's avatar
toy committed
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
	    (format t "p0, q0 = ~A ~A~%" p0 q0)
	    (format t "p1, q1 = ~A ~A~%" p1 q1)
	    (format t "p2, q2 = ~A ~A~%" p2 q2)
	    (format t "p3, q3 = ~A ~A~%" p3 q3)
	    (format t "p4, q4 = ~A ~A~%" p4 q4)
	    (multiple-value-setq (p1 p2 q0)
	      (three-sum p1 p2 q0))
	    (format t "p1 = ~A~%" p1)
	    (format t "p2 = ~A~%" p2)
	    (format t "q0 = ~A~%" q0)
	    ;; five-three-sum
	    (multiple-value-setq (p2 p3 p4)
	      (three-sum p2 p3 p4))
	    (format t "p2 = ~A~%" p2)
	    (format t "p3 = ~A~%" p3)
	    (format t "p4 = ~A~%" p4)
	    (multiple-value-setq (q1 q2)
632
	      (two-sum q1 q2))
toy's avatar
toy committed
633
	    (multiple-value-bind (s0 t0)
634
		(two-sum p2 q1)
toy's avatar
toy committed
635
	      (multiple-value-bind (s1 t1)
636
		  (two-sum p3 q2)
toy's avatar
toy committed
637
		(multiple-value-setq (s1 t0)
638
		  (two-sum s1 t0))
639
		(let ((s2 (cl:+ t0 t1 p4))
toy's avatar
toy committed
640
		      (p2 s0)
641
		      (p3 (cl:+ (cl:* (qd-2 a)
toy's avatar
toy committed
642
				(kernel:double-double-hi b))
643
			     (cl:* (qd-3 a)
toy's avatar
toy committed
644
645
				(kernel:double-double-lo b))
			     q3 q4)))
646
		  (multiple-value-setq (p3 q0)
toy's avatar
toy committed
647
		    (three-sum2 p3 q0 s1))
648
		  (let ((p4 (cl:+ q0 s2)))
toy's avatar
toy committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
		    (multiple-value-call #'%make-qd-d
		      (renorm-5 p0 p1 p2 p3 p4))))))))))))

;; a0 * b0                    0
;;      a0 * b1               1
;;      a1 * b0               2
;;           a0 * b2          3
;;           a1 * b1          4
;;           a2 * b0          5
;;                a0 * b3     6
;;                a1 * b2     7
;;                a2 * b1     8
;;                a3 * b0     9 

;; Works
;; (mul-qd (sqrt-qd (make-qd-dd 2w0 0w0)) (make-qd-dd 10w0 0w0)) ->
;; 14.1421356237309504880168872420969807856967187537694807317667973799q0
;;
;; Clisp says
;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
669
670

(defun mul-qd-t (a b target)
671
  (declare (type %quad-double a b #+oct-array target)
672
	   (optimize (speed 3)
Raymond Toy's avatar
Raymond Toy committed
673
		     (space 0))
674
	   (inline float-infinity-p)
675
	   #+(and cmu (not oct-array))
676
	   (ignore target))
Raymond Toy's avatar
Raymond Toy committed
677
678
679
680
681
682
  (with-qd-parts (a0 a1 a2 a3)
      a
    (declare (double-float a0 a1 a2 a3))
    (with-qd-parts (b0 b1 b2 b3)
	b
      (declare (double-float b0 b1 b2 b3))
toy's avatar
toy committed
683
      (multiple-value-bind (p0 q0)
684
	  (two-prod a0 b0)
685
686
	#+cmu
	(when (float-infinity-p p0)
687
	  (return-from mul-qd-t (%store-qd-d target p0 0d0 0d0 0d0)))
toy's avatar
toy committed
688
	(multiple-value-bind (p1 q1)
689
	    (two-prod a0 b1)
toy's avatar
toy committed
690
	  (multiple-value-bind (p2 q2)
691
	      (two-prod a1 b0)
toy's avatar
toy committed
692
	    (multiple-value-bind (p3 q3)
693
		(two-prod a0 b2)
toy's avatar
toy committed
694
	      (multiple-value-bind (p4 q4)
695
		  (two-prod a1 b1)
toy's avatar
toy committed
696
		(multiple-value-bind (p5 q5)
697
		    (two-prod a2 b0)
toy's avatar
toy committed
698
		  ;; Start accumulation
699
		  (three-sum p1 p2 q0 p1 p2 q0)
toy's avatar
toy committed
700
701

		  ;; six-three-sum of p2, q1, q2, p3, p4, p5
702
703
		  (three-sum p2 q1 q2 p2 q1 q2)
		  (three-sum p3 p4 p5 p3 p4 p5)
toy's avatar
toy committed
704
		  ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5)
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
		  (let* ((s0 0d0)
			 (s1 s0)
			 (t0 s0)
			 (t1 s0))
		    (declare (double-float s0 s1 t0 t1))
		    (two-sum s0 t0 p2 p3)
		    (two-sum s1 t1 q1 p4)
		    (let ((s2 (cl:+ q2 p5)))
		      (declare (double-float s2))
		      (two-sum s1 t0 s1 t0)
		      (incf s2 (cl:+ t0 t1))
		      ;; O(eps^3) order terms.  This is the sloppy
		      ;; multiplication version.  Should we use
		      ;; the precise version?  It's significantly
		      ;; more complex.
toy's avatar
toy committed
720
			  
721
722
723
724
725
726
727
728
729
730
731
		      (incf s1 (cl:+ (cl:* a0 b3)
				     (cl:* a1 b2)
				     (cl:* a2 b1)
				     (cl:* a3 b0)
				     q0 q3 q4 q5))
		      #+nil
		      (format t "p0,p1,s0,s1,s2 = ~a ~a ~a ~a ~a~%"
			      p0 p1 s0 s1 s2)
		      (multiple-value-bind (r0 r1 s0 s1)
			  (renorm-5 p0 p1 s0 s1 s2)
			(if (zerop r0)
732
733
734
			    (%store-qd-d target p0 0d0 0d0 0d0)
			    (%store-qd-d target r0 r1 s0 s1))))))))))))))

735
736
737
738
(defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
  "Returns the product of the %QUAD-DOUBLE numbers A and B.
If TARGET is given, TARGET is destructively modified to contain the result."
  (mul-qd-t a b target))
toy's avatar
toy committed
739
740
741
742
743
744
745
746
747
748
749
750
751

;; This is the non-sloppy version.  I think this works just fine, but
;; since qd defaults to the sloppy multiplication version, we do the
;; same.
#+nil
(defun mul-qd (a b)
  (declare (type %quad-double a b)
	   (optimize (speed 3)))
  (multiple-value-bind (a0 a1 a2 a3)
      (qd-parts a)
    (multiple-value-bind (b0 b1 b2 b3)
	(qd-parts b)
      (multiple-value-bind (p0 q0)
752
	  (two-prod a0 b0)
toy's avatar
toy committed
753
	(multiple-value-bind (p1 q1)
754
	    (two-prod a0 b1)
toy's avatar
toy committed
755
	  (multiple-value-bind (p2 q2)
756
	      (two-prod a1 b0)
toy's avatar
toy committed
757
	    (multiple-value-bind (p3 q3)
758
		(two-prod a0 b2)
toy's avatar
toy committed
759
	      (multiple-value-bind (p4 q4)
760
		  (two-prod a1 b1)
toy's avatar
toy committed
761
762
763
764
765
766
767
768
769
		(declare (double-float q4))
		#+nil
		(progn
		  (format t"p0, q0 = ~a ~a~%" p0 q0)
		  (format t"p1, q1 = ~a ~a~%" p1 q1)
		  (format t"p2, q2 = ~a ~a~%" p2 q2)
		  (format t"p3, q3 = ~a ~a~%" p3 q3)
		  (format t"p4, q4 = ~a ~a~%" p4 q4))
		(multiple-value-bind (p5 q5)
770
		    (two-prod a2 b0)
toy's avatar
toy committed
771
772
773
774
775
776
777
778
779
780
781
782
783
784
		  #+nil
		  (format t"p5, q5 = ~a ~a~%" p5 q5)
		  ;; Start accumulation
		  (multiple-value-setq (p1 p2 q0)
		    (three-sum p1 p2 q0))
		  #+nil
		  (format t "p1 p2 q0 = ~a ~a ~a~%" p1 p2 q0)
		  ;; six-three-sum of p2, q1, q2, p3, p4, p5
		  (multiple-value-setq (p2 q1 q2)
		    (three-sum p2 q1 q2))
		  (multiple-value-setq (p3 p4 p5)
		    (three-sum p3 p4 p5))
		  ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5)
		  (multiple-value-bind (s0 t0)
785
		      (two-sum p2 p3)
toy's avatar
toy committed
786
		    (multiple-value-bind (s1 t1)
787
			(two-sum q1 p4)
toy's avatar
toy committed
788
		      (declare (double-float t1))
789
		      (let ((s2 (cl:+ q2 p5)))
toy's avatar
toy committed
790
791
			(declare (double-float s2))
			(multiple-value-bind (s1 t0)
792
			    (two-sum s1 t0)
toy's avatar
toy committed
793
			  (declare (double-float s1))
794
			  (incf s2 (cl:+ t0 t1))
toy's avatar
toy committed
795
			  (multiple-value-bind (p6 q6)
796
			      (two-prod a0 b3)
toy's avatar
toy committed
797
			    (multiple-value-bind (p7 q7)
798
				(two-prod a1 b2)
toy's avatar
toy committed
799
			      (multiple-value-bind (p8 q8)
800
				  (two-prod a2 b1)
toy's avatar
toy committed
801
				(multiple-value-bind (p9 q9)
802
				    (two-prod a3 b0)
toy's avatar
toy committed
803
				  (multiple-value-setq (q0 q3)
804
				    (two-sum q0 q3))
toy's avatar
toy committed
805
				  (multiple-value-setq (q4 q5)
806
				    (two-sum q4 q5))
toy's avatar
toy committed
807
				  (multiple-value-setq (p6 p7)
808
				    (two-sum p6 p7))
toy's avatar
toy committed
809
				  (multiple-value-setq (p8 p9)
810
				    (two-sum p8 p9))
toy's avatar
toy committed
811
812
				  ;; (t0,t1) = (q0,q3)+(q4,q5)
				  (multiple-value-setq (t0 t1)
813
				    (two-sum q0 q4))
814
				  (setf t1 (cl:+ q3 q5))
toy's avatar
toy committed
815
816
				  ;; (r0,r1) = (p6,p7)+(p8,p9)
				  (multiple-value-bind (r0 r1)
817
				      (two-sum p6 p8)
toy's avatar
toy committed
818
				    (declare (double-float r1))
819
				    (incf r1 (cl:+ p7 p9))
toy's avatar
toy committed
820
821
				    ;; (q3,q4) = (t0,t1)+(r0,r1)
				    (multiple-value-setq (q3 q4)
822
				      (two-sum t0 r0))
823
				    (incf q4 (cl:+ t1 r1))
toy's avatar
toy committed
824
825
				    ;; (t0,t1)=(q3,q4)+s1
				    (multiple-value-setq (t0 t1)
826
				      (two-sum q3 s1))
toy's avatar
toy committed
827
828
829
				    (incf t1 q4)
				    ;; O(eps^4) terms
				    (incf t1
830
831
832
					  (cl:+ (cl:* a1 b3)
					     (cl:* a2 b2)
					     (cl:* a3 b1)
toy's avatar
toy committed
833
834
835
836
837
838
839
840
					     q6 q7 q8 q9
					     s2))
				    #+nil
				    (format t "p0,p1,s0,t0,t1 = ~a ~a ~a ~a ~a~%"
					    p0 p1 s0 t0 t1)
				    (multiple-value-call #'%make-qd-d
				      (renorm-5 p0 p1 s0 t0 t1))))))))))))))))))))

841
(defun sqr-qd-t (a target)
toy's avatar
toy committed
842
  "Square A"
843
  (declare (type %quad-double a #+oct-array target)
844
	   (optimize (speed 3)
845
		     (space 0))
846
	   #+(and cmu (not oct-array))
847
	   (ignore target))
toy's avatar
toy committed
848
  (multiple-value-bind (p0 q0)
849
      (two-sqr (qd-0 a))
toy's avatar
toy committed
850
    (multiple-value-bind (p1 q1)
851
	(two-prod (cl:* 2 (qd-0 a)) (qd-1 a))
toy's avatar
toy committed
852
      (multiple-value-bind (p2 q2)
853
	  (two-prod (cl:* 2 (qd-0 a)) (qd-2 a))
toy's avatar
toy committed
854
	(multiple-value-bind (p3 q3)
855
	    (two-sqr (qd-1 a))
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
	  (two-sum p1 q0 q0 p1)
	  (two-sum q0 q1 q0 q1)
	  (two-sum p2 p3 p2 p3)

	  (let* ((s0 0d0)
		 (t0 s0)
		 (s1 s0)
		 (t1 s0))
	    (declare (double-float s0 s1 t0 t1))
	    (two-sum s0 t0 q0 p2)
	    (two-sum s1 t1 q1 p3)
	    (two-sum s1 t0 s1 t0)
	    (incf t0 t1)
	    (quick-two-sum s1 t0 s1 t0)
	    (quick-two-sum p2 t1 s0 s1)
	    (quick-two-sum p3 q0 t1 t0)

	    (let ((p4 (cl:* 2 (qd-0 a) (qd-3 a)))
		  (p5 (cl:* 2 (qd-1 a) (qd-2 a))))
	      (declare (double-float p4 p5))
	      (two-sum p4 p5 p4 p5)
	      (two-sum q2 q3 q2 q3)
	      (two-sum t0 t1 p4 q2)

	      (incf t1 (cl:+ p5 q3))

	      (two-sum p3 p4 p3 t0)
	      (incf p4 (cl:+ q0 t1))

885
886
	      (multiple-value-bind (a0 a1 a2 a3)
		  (renorm-5 p0 p1 p2 p3 p4)
887
		(%store-qd-d target a0 a1 a2 a3)))))))))
toy's avatar
toy committed
888

889
890
891
892
(defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
  "Return the square of the %QUAD-DOUBLE number A.  If TARGET is also given,
it is destructively modified with the result."
  (sqr-qd-t a target))
893
894
895

#+nil
(defun div-qd-t (a b target)
896
897
  (declare (type %quad-double a b)
	   (optimize (speed 3)
Raymond Toy's avatar
Raymond Toy committed
898
		     (space 0))
899
900
901
	   (inline float-infinity-p)
	   #+cmu
	   (ignore target))
toy's avatar
toy committed
902
903
  (let ((a0 (qd-0 a))
	(b0 (qd-0 b)))
904
    (let* ((q0 (cl:/ a0 b0))
toy's avatar
toy committed
905
	   (r (sub-qd a (mul-qd-d b q0)))
906
	   (q1 (cl:/ (qd-0 r) b0)))
907

908
      (when (float-infinity-p q0)
909
	(return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0)))
toy's avatar
toy committed
910
      (setf r (sub-qd r (mul-qd-d b q1)))
911
      (let ((q2 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
912
	(setf r (sub-qd r (mul-qd-d b q2)))
913
	(let ((q3 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
914
915
	  (multiple-value-bind (q0 q1 q2 q3)
	      (renorm-4 q0 q1 q2 q3)
916
917
918
	    (%store-qd-d target q0 q1 q2 q3)))))))

(defun div-qd-t (a b target)
919
  (declare (type %quad-double a b #+oct-array target)
920
921
922
	   (optimize (speed 3)
		     (space 0))
	   (inline float-infinity-p)
923
	   #+(and cmu (not oct-array))
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
	   (ignore target))
  (let ((a0 (qd-0 a))
	(b0 (qd-0 b)))
    (let* ((q0 (cl:/ a0 b0))
	   (r (%make-qd-d 0d0 0d0 0d0 0d0)))
      (mul-qd-d b q0 r)
      (sub-qd a r r)
      (let* ((q1 (cl:/ (qd-0 r) b0))
	     (temp (mul-qd-d b q1)))
	(when (float-infinity-p q0)
	  (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0)))
	
	(sub-qd r temp r)
	(let ((q2 (cl:/ (qd-0 r) b0)))
	  (mul-qd-d b q2 temp)
	  (sub-qd r temp r)
	  (let ((q3 (cl:/ (qd-0 r) b0)))
	    (multiple-value-bind (q0 q1 q2 q3)
		(renorm-4 q0 q1 q2 q3)
	      (%store-qd-d target q0 q1 q2 q3))))))))
toy's avatar
toy committed
944

945
946
947
948
949
(defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
  "Return the quotient of the two %QUAD-DOUBLE numbers A and B.
If TARGET is given, it destrutively modified with the result."
  (div-qd-t a b target))

950
951
(declaim (inline invert-qd))

952
953
(defun invert-qd (v)
  ;; a quartic newton iteration for 1/v
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
  ;; to invert v, start with a good guess, x.
  ;; let h= 1-v*x  ;; h is small
  ;; return x+ x*(h+h^2+h^3) . compute h3 in double-float
  ;; enough accuracy.
    
  (let* 
      ((x (%make-qd-d (cl:/ (qd-0 v)) 0d0 0d0 0d0))
       (h (add-qd-d (neg-qd (mul-qd v x))
		    1.0d0))
       (h2 (mul-qd h h)) ;also use h2 for target
       (h3 (* (qd-0 h) (qd-0 h2))))
    (declare (type %quad-double v h h2)
	     (double-float h3))
    (add-qd x
	    (mul-qd x
		    (add-qd h (add-qd-d h2 h3))))))

(defun div-qd-i (a b)
  (declare (type %quad-double a b)
	   (optimize (speed 3)
		     (space 0)))
  (mul-qd a (invert-qd b)))
  
toy's avatar
toy committed
977
978
979
980
981
982
;; Non-sloppy divide
#+(or)
(defun div-qd (a b)
  (declare (type %quad-double a b))
  (let ((a0 (qd-0 a))
	(b0 (qd-0 b)))
983
    (let* ((q0 (cl:/ a0 b0))
toy's avatar
toy committed
984
	   (r (sub-qd a (mul-qd-d b q0)))
985
	   (q1 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
986
      (setf r (sub-qd r (mul-qd-d b q1)))
987
      (let ((q2 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
988
	(setf r (sub-qd r (mul-qd-d b q2)))
989
	(let ((q3 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
990
	  (setf r (sub-qd r (mul-qd-d b q3)))
991
	  (let ((q4 (cl:/ (qd-0 r) b0)))
toy's avatar
toy committed
992
993
994
995
996
997
	  (multiple-value-bind (q0 q1 q2 q3)
	      (renorm-5 q0 q1 q2 q3 q4)
	    (%make-qd-d q0 q1 q2 q3))))))))

;; quad-double / double
(defun div-qd-d (a b)
Raymond Toy's avatar
Raymond Toy committed
998
  "Divides the %QUAD-DOUBLE number A by the DOUBLE-FLOAT number B."
toy's avatar
toy committed
999
1000
  (declare (type %quad-double a)
	   (double-float b)
1001
	   (optimize (speed 3)
Raymond Toy's avatar
Raymond Toy committed
1002
		     (space 0))
1003
	   (inline float-infinity-p))
toy's avatar
toy committed
1004
1005
1006
  ;; Compute approximate quotient using high order doubles, then
  ;; correct it 3 times using the remainder.  Analogous to long
  ;; division.
1007
  (let ((q0 (cl:/ (qd-0 a) b)))
1008

1009
1010
1011
    (when (float-infinity-p q0)
      (return-from div-qd-d (%make-qd-d q0 0d0 0d0 0d0)))

toy's avatar
toy committed
1012
1013
    ;; Compute remainder a - q0 * b
    (multiple-value-bind (t0 t1)
1014
	(two-prod q0 b)
1015
1016
      (let ((r #+cmu (sub-qd-dd a (kernel:make-double-double-float t0 t1))
	       #-cmu (sub-qd a (make-qd-d t0 t1 0d0 0d0))))
toy's avatar
toy committed
1017
	;; First correction
1018
	(let ((q1 (cl:/ (qd-0 r) b)))
toy's avatar
toy committed
1019
	  (multiple-value-bind (t0 t1)