Context Navigation

root/qd.lisp @ 258cd829289826458d6075caa279b253b8848749

Revision 258cd829289826458d6075caa279b253b8848749, 39.8 KB (checked in by Raymond Toy <toy.raymond@…>, 2 years ago)
Remove extra blank lines. Fix #1.
Property mode set to `100644`

Line
1	;;;; -- Mode: lisp --
2	;;;;
3	;;;; Copyright (c) 2007, 2008, 2011 Raymond Toy
4	;;;;
5	;;;; Permission is hereby granted, free of charge, to any person
6	;;;; obtaining a copy of this software and associated documentation
7	;;;; files (the "Software"), to deal in the Software without
8	;;;; restriction, including without limitation the rights to use,
9	;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10	;;;; copies of the Software, and to permit persons to whom the
11	;;;; Software is furnished to do so, subject to the following
12	;;;; conditions:
13	;;;;
14	;;;; The above copyright notice and this permission notice shall be
15	;;;; included in all copies or substantial portions of the Software.
16	;;;;
17	;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18	;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19	;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20	;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21	;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22	;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23	;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24	;;;; OTHER DEALINGS IN THE SOFTWARE.
25
26	;;; This file contains the core routines for basic arithmetic
27	;;; operations on a %quad-double. This includes addition,
28	;;; subtraction, multiplication, division, negation. and absolute
29	;;; value. Some additional helper functions are included such as
30	;;; raising to an integer power. and the n'th root of a (non-negative)
31	;;; %quad-double. The basic comparison operators are included, and
32	;;; some simple tests for zerop, onep, plusp, and minusp.
33	;;;
34	;;; The basic algorithms are based on Yozo Hida's double-double
35	;;; implementation. However, some were copied from CMUCL and modified
36	;;; to support quad-doubles.
37
38
39	(in-package #:octi)
40
41	#+cmu
42	(eval-when (:compile-toplevel)
43	(setf ext:inline-expansion-limit 1600))
44
45	;; All of the following functions should be inline.
46	;;(declaim (inline three-sum2))
47
48	;; Internal routines for implementing quad-double.
49
50	#+cmu
51	(defmacro two-sum (s e x y)
52	`(multiple-value-setq (,s ,e)
53	(c::two-sum ,x ,y)))
54
55
56	(defmacro three-sum (s0 s1 s2 a b c)
57	(let ((aa (gensym))
58	(bb (gensym))
59	(cc (gensym))
60	(t1 (gensym))
61	(t2 (gensym))
62	(t3 (gensym)))
63	`(let* ((,aa ,a)
64	(,bb ,b)
65	(,cc ,c)
66	(,t1 0d0)
67	(,t2 ,t1)
68	(,t3 ,t1))
69	(declare (double-float ,t1 ,t2 ,t3))
70	(two-sum ,t1 ,t2 ,aa ,bb)
71	(two-sum ,s0 ,t3 ,cc ,t1)
72	(two-sum ,s1 ,s2 ,t2 ,t3))))
73
74
75
76	#+nil
77	(defun three-sum2 (a b c)
78	(declare (double-float a b c)
79	(optimize (speed 3)))
80	(let* ((t1 0d0)
81	(t2 t1)
82	(t3 t1))
83	(two-sum t1 t2 a b)
84	(two-sum a t3 c t1)
85	(values a (cl:+ t2 t3))))
86
87	(defmacro three-sum2 (s0 s1 a b c)
88	(let ((aa (gensym))
89	(bb (gensym))
90	(cc (gensym))
91	(t1 (gensym))
92	(t2 (gensym))
93	(t3 (gensym)))
94	`(let* ((,aa ,a)
95	(,bb ,b)
96	(,cc ,c)
97	(,t1 0d0)
98	(,t2 ,t1)
99	(,t3 ,t1))
100	(declare (double-float ,t1 ,t2 ,t3))
101	(two-sum ,t1 ,t2 ,aa ,bb)
102	(two-sum ,s0 ,t3 ,cc ,t1)
103	(setf ,s1 (+ ,t2 ,t3)))))
104
105	;; Not needed????
106	#+nil
107	(declaim (inline quick-three-accum))
108	#+nil
109	(defun quick-three-accum (a b c)
110	(declare (double-float a b c)
111	(optimize (speed 3) (space 0)))
112	(multiple-value-bind (s b)
113	(two-sum b c)
114	(multiple-value-bind (s a)
115	(two-sum a s)
116	(let ((za (/= a 0))
117	(zb (/= b 0)))
118	(when (and za zb)
119	(return-from quick-three-accum (values (cl:+ s 0d0) (cl:+ a 0d0) (cl:+ b 0d0))))
120	(when (not za)
121	(setf a s)
122	(setf s 0d0))
123	(when (not zb)
124	(setf b a)
125	(setf a s))
126	(values 0d0 a b)))))
127
128
129	;; These functions are quite short, so we inline them to minimize
130	;; consing.
131	(declaim (inline make-qd-d
132	add-d-qd
133	add-dd-qd
134	neg-qd
135	neg-qd-t
136	sub-qd
137	sub-qd-dd
138	sub-qd-d
139	sub-d-qd
140	#+cmu
141	make-qd-dd
142	abs-qd
143	qd->
144	qd-<
145	qd->=
146	qd-<=
147	zerop-qd
148	onep-qd
149	plusp-qd
150	minusp-qd
151	qd-=
152	scale-float-qd))
153
154	;; Should these functions be inline? The QD C++ source has these as
155	;; inline functions, but these are fairly large functions. However,
156	;; inlining makes quite a big difference in speed and consing.
157	#+cmu
158	(declaim (#+qd-inline inline
159	#-qd-inline ext:maybe-inline
160	renorm-4
161	renorm-5
162	add-qd-d
163	add-qd-d-t
164	add-qd-dd
165	add-qd
166	add-qd-t
167	mul-qd-d
168	mul-qd-d-t
169	mul-qd-dd
170	mul-qd
171	mul-qd-t
172	sqr-qd
173	sqr-qd-t
174	div-qd
175	div-qd-t
176	div-qd-d
177	div-qd-dd))
178
179	#+cmu
180	(defmacro quick-two-sum (s e x y)
181	`(multiple-value-setq (,s ,e)
182	(c::quick-two-sum ,x ,y)))
183
184	#-(or qd-inline (not cmu))
185	(declaim (ext:start-block renorm-4 renorm-5
186	make-qd-d
187	add-qd-d add-d-qd add-qd-dd
188	add-dd-qd
189	add-qd add-qd-t
190	add-qd-d-t
191	neg-qd
192	sub-qd sub-qd-dd sub-qd-d sub-d-qd
193	mul-qd-d mul-qd-dd
194	mul-qd
195	mul-qd-t
196	mul-qd-d-t
197	sqr-qd
198	sqr-qd-t
199	div-qd div-qd-d div-qd-dd
200	div-qd-t
201	make-qd-dd
202	))
203
204	#+(or)
205	(defun quick-renorm (c0 c1 c2 c3 c4)
206	(declare (double-float c0 c1 c2 c3 c4)
207	(optimize (speed 3)))
208	(multiple-value-bind (s t3)
209	(quick-two-sum c3 c4)
210	(multiple-value-bind (s t2)
211	(quick-two-sum c2 s)
212	(multiple-value-bind (s t1)
213	(quick-two-sum c1 s)
214	(multiple-value-bind (c0 t0)
215	(quick-two-sum c0 s)
216	(multiple-value-bind (s t2)
217	(quick-two-sum t2 t3)
218	(multiple-value-bind (s t1)
219	(quick-two-sum t1 s)
220	(multiple-value-bind (c1 t0)
221	(quick-two-sum t0 s)
222	(multiple-value-bind (s t1)
223	(quick-two-sum t1 t2)
224	(multiple-value-bind (c2 t0)
225	(quick-two-sum t0 s)
226	(values c0 c1 c2 (cl:+ t0 t1))))))))))))
227
228	(defun renorm-4 (c0 c1 c2 c3)
229	(declare (double-float c0 c1 c2 c3)
230	(optimize (speed 3) (safety #-allegro 0 #+allegro 1) (debug 0)))
231	(let ((s0 0d0)
232	(s1 0d0)
233	(s2 0d0)
234	(s3 0d0))
235	(declare (double-float s0 s1 s2 s3))
236	(quick-two-sum s0 c3 c2 c3)
237	(quick-two-sum s0 c2 c1 s0)
238	(quick-two-sum c0 c1 c0 s0)
239	(setf s0 c0)
240	(setf s1 c1)
241	(cond ((/= s1 0)
242	(quick-two-sum s1 s2 s1 c2)
243	(if (/= s2 0)
244	(quick-two-sum s2 s3 s2 c3)
245	(quick-two-sum s1 s2 s1 c3)))
246	(t
247	(quick-two-sum s0 s1 s0 c2)
248	(if (/= s1 0)
249	(quick-two-sum s1 s2 s1 c3)
250	(quick-two-sum s0 s1 s0 c3))))
251	(values s0 s1 s2 s3)))
252
253	(defun renorm-5 (c0 c1 c2 c3 c4)
254	(declare (double-float c0 c1 c2 c3 c4)
255	(optimize (speed 3) (safety #-allegro 0 #+allegro 1)))
256	(let ((s0 0d0)
257	(s1 0d0)
258	(s2 0d0)
259	(s3 0d0))
260	(declare (double-float s0 s1 s2 s3))
261	(quick-two-sum s0 c4 c3 c4)
262	(quick-two-sum s0 c3 c2 s0)
263	(quick-two-sum s0 c2 c1 s0)
264	(quick-two-sum c0 c1 c0 s0)
265	(quick-two-sum s0 s1 c0 c1)
266	(cond ((/= s1 0)
267	(quick-two-sum s1 s2 s1 c2)
268	(cond ((/= s2 0)
269	(quick-two-sum s2 s3 s2 c3)
270	(if (/= s3 0)
271	(incf s3 c4)
272	(incf s2 c4)))
273	(t
274	(quick-two-sum s1 s2 s1 c3)
275	(if (/= s2 0)
276	(quick-two-sum s2 s3 s2 c4)
277	(quick-two-sum s1 s2 s1 c4)))))
278	(t
279	(quick-two-sum s0 s1 s0 c2)
280	(cond ((/= s1 0)
281	(quick-two-sum s1 s2 s1 c3)
282	(if (/= s2 0)
283	(quick-two-sum s2 s3 s2 c4)
284	(quick-two-sum s1 s2 s1 c4)))
285	(t
286	(quick-two-sum s0 s1 s0 c3)
287	(if (/= s1 0)
288	(quick-two-sum s1 s2 s1 c4)
289	(quick-two-sum s0 s1 s0 c4))))))
290	(values s0 s1 s2 s3)))
291
292	(defun make-qd-d (a0 &optional (a1 0d0 a1-p) (a2 0d0) (a3 0d0))
293	"Create a %QUAD-DOUBLE from four double-floats, appropriately
294	normalizing the result from the four double-floats. A0 is the most
295	significant part of the %QUAD-DOUBLE, and A3 is the least. The optional
296	parameters default to 0.
297	"
298	(declare (double-float a0 a1 a2 a3)
299	(optimize (speed 3)
300	#+cmu
301	(ext:inhibit-warnings 3)))
302	(if a1-p
303	(multiple-value-bind (s0 s1 s2 s3)
304	(renorm-4 a0 a1 a2 a3)
305	(%make-qd-d s0 s1 s2 s3))
306	(%make-qd-d a0 0d0 0d0 0d0)))
307
308	;;;; Addition
309
310	;; Quad-double + double
311	(defun add-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
312	"Return the sum of the %QUAD-DOUBLE A and the DOUBLE-FLOAT B.
313	If TARGET is given, TARGET is destructively modified to contain the result."
314	(add-qd-d-t a b target))
315
316	(defun add-qd-d-t (a b target)
317	"Add a quad-double A and a double-float B"
318	(declare (type %quad-double a #+oct-array target)
319	(double-float b)
320	(optimize (speed 3)
321	(space 0))
322	(inline float-infinity-p)
323	#+(and cmu (not oct-array)) (ignore target))
324	(let* ((c0 0d0)
325	(e c0)
326	(c1 c0)
327	(c2 c0)
328	(c3 c0))
329	(declare (double-float e c0 c1 c2 c3))
330	(two-sum c0 e (qd-0 a) b)
331
332	(when (float-infinity-p c0)
333	(return-from add-qd-d-t (%store-qd-d target c0 0d0 0d0 0d0)))
334	(two-sum c1 e (qd-1 a) e)
335	(two-sum c2 e (qd-2 a) e)
336	(two-sum c3 e (qd-3 a) e)
337	(multiple-value-bind (r0 r1 r2 r3)
338	(renorm-5 c0 c1 c2 c3 e)
339	(if (and (zerop (qd-0 a)) (zerop b))
340	(%store-qd-d target c0 0d0 0d0 0d0)
341	(%store-qd-d target r0 r1 r2 r3)))))
342
343	(defun add-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
344	"Return the sum of the DOUBLE-FLOAT A and the %QUAD-DOUBLE B.
345	If TARGET is given, TARGET is destructively modified to contain the result."
346	(declare (double-float a)
347	(type %quad-double b)
348	(optimize (speed 3))
349	#+(and cmu (not oct-array)) (ignore target))
350	(add-qd-d b a #+oct-array target))
351
352	#+cmu
353	(defun add-qd-dd (a b)
354	"Add a quad-double A and a double-double B"
355	(declare (type %quad-double a)
356	(double-double-float b)
357	(optimize (speed 3)
358	(space 0)))
359	(let* ((s0 0d0)
360	(t0 s0)
361	(s1 s0)
362	(t1 s0)
363	(s2 s0)
364	(s3 s0))
365	(declare (double-float s0 s1 s3 t0 t1 s2))
366	(two-sum s0 t1 (qd-0 a) (kernel:double-double-hi b))
367	(two-sum s1 t1 (qd-1 a) (kernel:double-double-lo b))
368	(two-sum s1 t0 s1 t0)
369	(three-sum s2 t0 t1 (qd-2 a) t0 t1)
370	(two-sum s3 t0 t0 (qd-3 a))
371	(incf t0 t1)
372	(multiple-value-bind (a0 a1 a2 a3)
373	(renorm-5 s0 s1 s2 s3 t0)
374	(declare (double-float a0 a1 a2 a3))
375	(%make-qd-d a0 a1 a2 a3))))
376
377	#+cmu
378	(defun add-dd-qd (a b)
379	(declare (double-double-float a)
380	(type %quad-double b)
381	(optimize (speed 3)
382	(space 0)))
383	(add-qd-dd b a))
384
385
386	#+(or)
387	(defun add-qd-1 (a b)
388	(declare (type %quad-double a b)
389	(optimize (speed 3)))
390	(multiple-value-bind (s0 t0)
391	(two-sum (qd-0 a) (qd-0 b))
392	(multiple-value-bind (s1 t1)
393	(two-sum (qd-1 a) (qd-1 b))
394	(multiple-value-bind (s2 t2)
395	(two-sum (qd-2 a) (qd-2 b))
396	(multiple-value-bind (s3 t3)
397	(two-sum (qd-3 a) (qd-3 b))
398	(multiple-value-bind (s1 t0)
399	(two-sum s1 t0)
400	(multiple-value-bind (s2 t0 t1)
401	(three-sum s2 t0 t1)
402	(multiple-value-bind (s3 t0)
403	(three-sum2 s3 t0 t2)
404	(let ((t0 (cl:+ t0 t1 t3)))
405	(multiple-value-call #'%make-qd-d
406	(renorm-5 s0 s1 s2 s3 t0)))))))))))
407
408	;; Same as above, except that everything is expanded out for compilers
409	;; which don't do a very good job with dataflow. CMUCL is one of
410	;; those compilers.
411
412	(defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
413	"Return the sum of the %QUAD-DOUBLE numbers A and B.
414	If TARGET is given, TARGET is destructively modified to contain the result."
415	(add-qd-t a b target))
416
417
418	(defun add-qd-t (a b target)
419	(declare (type %quad-double a b #+oct-array target)
420	(optimize (speed 3)
421	(space 0))
422	#+(and cmu (not oct-array))
423	(ignore target))
424	;; This is the version that is NOT IEEE. Should we use the IEEE
425	;; version? It's quite a bit more complicated.
426	;;
427	;; In addition, this is reorganized to minimize data dependency.
428	(with-qd-parts (a0 a1 a2 a3)
429	a
430	(declare (double-float a0 a1 a2 a3))
431	(with-qd-parts (b0 b1 b2 b3)
432	b
433	(declare (double-float b0 b1 b2 b3))
434	(let ((s0 (cl:+ a0 b0))
435	(s1 (cl:+ a1 b1))
436	(s2 (cl:+ a2 b2))
437	(s3 (cl:+ a3 b3)))
438	(declare (double-float s0 s1 s2 s3)
439	(inline float-infinity-p))
440
441	(when (float-infinity-p s0)
442	(return-from add-qd-t (%store-qd-d target s0 0d0 0d0 0d0)))
443	(let ((v0 (cl:- s0 a0))
444	(v1 (cl:- s1 a1))
445	(v2 (cl:- s2 a2))
446	(v3 (cl:- s3 a3)))
447	(declare (double-float v0 v1 v2 v3))
448	(let ((u0 (cl:- s0 v0))
449	(u1 (cl:- s1 v1))
450	(u2 (cl:- s2 v2))
451	(u3 (cl:- s3 v3)))
452	(declare (double-float u0 u1 u2 u3))
453	(let ((w0 (cl:- a0 u0))
454	(w1 (cl:- a1 u1))
455	(w2 (cl:- a2 u2))
456	(w3 (cl:- a3 u3)))
457	(declare (double-float w0 w1 w2 w3))
458	(let ((u0 (cl:- b0 v0))
459	(u1 (cl:- b1 v1))
460	(u2 (cl:- b2 v2))
461	(u3 (cl:- b3 v3)))
462	(declare (double-float u0 u1 u2 u3))
463	(let ((t0 (cl:+ w0 u0))
464	(t1 (cl:+ w1 u1))
465	(t2 (cl:+ w2 u2))
466	(t3 (cl:+ w3 u3)))
467	(declare (double-float t0 t1 t2 t3))
468	(two-sum s1 t0 s1 t0)
469	(three-sum s2 t0 t1 s2 t0 t1)
470	(three-sum2 s3 t0 s3 t0 t2)
471	(setf t0 (cl:+ t0 t1 t3))
472	;; Renormalize
473	(multiple-value-setq (s0 s1 s2 s3)
474	(renorm-5 s0 s1 s2 s3 t0))
475	(if (and (zerop a0) (zerop b0))
476	(%store-qd-d target (+ a0 b0) 0d0 0d0 0d0)
477	(%store-qd-d target s0 s1 s2 s3)))))))))))
478
479	(defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
480	"Return the negative of the %QUAD-DOUBLE number A.
481	If 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.
494	If 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
509	the DOUBLE-FLOAT number B.
510
511	If 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
519	the %QUAD-DOUBLE number B.
520
521	If 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
538	the DOUBLE-FLOAT number B.
539
540	If 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.
669	If 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,
841	it 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.
895	If 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,
1054	which 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
1076	is 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 browser.

Download in other formats:

Original Format