Context Navigation

root/qd-methods.lisp @ 1e343ca0be49b9e56938a8063b4d2e679e338965

Revision 1e343ca0be49b9e56938a8063b4d2e679e338965, 32.0 KB (checked in by Raymond Toy <toy.raymond@…>, 2 years ago)

Fix issue with float-sign calling qfloat-sign.

If the optional arg to FLOAT-SIGN is not given, don't call qfloat-sign
with a second arg of NIL. This breaks things. Call qfloat-sign with
the same number of arguments as float-sign.

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	(in-package #:oct)
27
28	(defconstant +pi+
29	(make-instance 'qd-real :value octi:+qd-pi+)
30	"Pi represented as a QD-REAL")
31
32	(defconstant +pi/2+
33	(make-instance 'qd-real :value octi:+qd-pi/2+)
34	"Pi/2 represented as a QD-REAL")
35
36	(defconstant +pi/4+
37	(make-instance 'qd-real :value octi:+qd-pi/4+)
38	"Pi/4 represented as a QD-REAL")
39
40	(defconstant +2pi+
41	(make-instance 'qd-real :value octi:+qd-2pi+)
42	"2*pi represented as a QD-REAL")
43
44	(defconstant +log2+
45	(make-instance 'qd-real :value octi:+qd-log2+)
46	"Natural log of 2 represented as a QD-REAL")
47
48	;; How do we represent infinity for a QD-REAL? For now, we just make
49	;; the QD-REAL whose most significant part is infinity. Currently
50	;; only supported on CMUCL.
51	#+cmu
52	(defconstant +quad-double-float-positive-infinity+
53	(make-instance 'qd-real :value (make-qd-d ext:double-float-positive-infinity))
54	"One representation of positive infinity for QD-REAL")
55
56	#+cmu
57	(defconstant +quad-double-float-negative-infinity+
58	(make-instance 'qd-real :value (make-qd-d ext:double-float-negative-infinity))
59	"One representation of negative infinity for QD-REAL")
60
61	(defconstant +most-positive-quad-double-float+
62	(make-instance 'qd-real
63	:value (octi::%make-qd-d most-positive-double-float
64	(cl:scale-float most-positive-double-float (cl:* 1 -53))
65	(cl:scale-float most-positive-double-float (cl:* 2 -53))
66	(cl:scale-float most-positive-double-float (cl:* 3 -53))))
67	"Most positive representable QD-REAL")
68
69	(defconstant +least-positive-quad-double-float+
70	(make-instance 'qd-real
71	:value (make-qd-d least-positive-double-float))
72	"Least positive QD-REAL")
73
74	;; Not sure this is 100% correct, but I think if the first component
75	;; is any smaller than this, the last component would no longer be a
76	;; normalized double-float.
77	(defconstant +least-positive-normalized-quad-double-float+
78	(make-instance 'qd-real
79	:value (make-qd-d (cl:scale-float least-positive-normalized-double-float (cl:* 3 53))))
80	"Least positive normalized QD-REAL")
81
82	(defconstant +qd-real-one+
83	(make-instance 'qd-real :value (make-qd-d 1d0))
84	"QD-REAL representation of 1")
85
86
87	(defmethod make-qd ((x cl:rational))
88	;; We should do something better than this.
89	(make-instance 'qd-real :value (rational-to-qd x)))
90
91
92	(defmethod add1 ((a number))
93	(cl::1+ a))
94
95	(defmethod add1 ((a qd-real))
96	(make-instance 'qd-real :value (add-qd-d (qd-value a) 1d0)))
97
98	(defmethod sub1 ((a number))
99	(cl::1- a))
100
101	(defmethod sub1 ((a qd-real))
102	(make-instance 'qd-real :value (sub-qd-d (qd-value a) 1d0)))
103
104	(declaim (inline 1+ 1-))
105
106	(defun 1+ (x)
107	(add1 x))
108
109	(defun 1- (x)
110	(sub1 x))
111
112	(defmethod two-arg-+ ((a qd-real) (b qd-real))
113	(make-instance 'qd-real :value (add-qd (qd-value a) (qd-value b))))
114
115	(defmethod two-arg-+ ((a qd-real) (b cl:float))
116	(make-instance 'qd-real :value (add-qd-d (qd-value a) (cl:float b 1d0))))
117
118	#+cmu
119	(defmethod two-arg-+ ((a qd-real) (b ext:double-double-float))
120	(make-instance 'qd-real :value (add-qd-dd (qd-value a) b)))
121
122	(defmethod two-arg-+ ((a real) (b qd-real))
123	(+ b a))
124
125	(defmethod two-arg-+ ((a number) (b number))
126	(cl:+ a b))
127
128	(defun + (&rest args)
129	(if (null args)
130	0
131	(do ((args (cdr args) (cdr args))
132	(res (car args)
133	(two-arg-+ res (car args))))
134	((null args) res))))
135
136	(defmethod two-arg-- ((a qd-real) (b qd-real))
137	(make-instance 'qd-real :value (sub-qd (qd-value a) (qd-value b))))
138
139	(defmethod two-arg-- ((a qd-real) (b cl:float))
140	(make-instance 'qd-real :value (sub-qd-d (qd-value a) (cl:float b 1d0))))
141
142	#+cmu
143	(defmethod two-arg-- ((a qd-real) (b ext:double-double-float))
144	(make-instance 'qd-real :value (sub-qd-dd (qd-value a) b)))
145
146	(defmethod two-arg-- ((a cl:float) (b qd-real))
147	(make-instance 'qd-real :value (sub-d-qd (cl:float a 1d0) (qd-value b))))
148
149	(defmethod two-arg-- ((a number) (b number))
150	(cl:- a b))
151
152	(defmethod unary-minus ((a number))
153	(cl:- a))
154
155	(defmethod unary-minus ((a qd-real))
156	(make-instance 'qd-real :value (neg-qd (qd-value a))))
157
158	(defun - (number &rest more-numbers)
159	(if more-numbers
160	(do ((nlist more-numbers (cdr nlist))
161	(result number))
162	((atom nlist) result)
163	(declare (list nlist))
164	(setq result (two-arg-- result (car nlist))))
165	(unary-minus number)))
166
167
168	(defmethod two-arg-* ((a qd-real) (b qd-real))
169	(make-instance 'qd-real :value (mul-qd (qd-value a) (qd-value b))))
170
171	(defmethod two-arg-* ((a qd-real) (b cl:float))
172	(make-instance 'qd-real :value (mul-qd-d (qd-value a) (cl:float b 1d0))))
173
174	#+cmu
175	(defmethod two-arg-* ((a qd-real) (b ext:double-double-float))
176	;; We'd normally want to use mul-qd-dd, but mul-qd-dd is broken.
177	(make-instance 'qd-real :value (mul-qd (qd-value a)
178	(make-qd-dd b 0w0))))
179
180	(defmethod two-arg-* ((a real) (b qd-real))
181	(* b a))
182
183	(defmethod two-arg-* ((a number) (b number))
184	(cl:* a b))
185
186	(defun * (&rest args)
187	(if (null args)
188	1
189	(do ((args (cdr args) (cdr args))
190	(res (car args)
191	(two-arg-* res (car args))))
192	((null args) res))))
193
194	(defmethod two-arg-/ ((a qd-real) (b qd-real))
195	(make-instance 'qd-real :value (div-qd (qd-value a) (qd-value b))))
196
197	(defmethod two-arg-/ ((a qd-real) (b cl:float))
198	(make-instance 'qd-real :value (div-qd-d (qd-value a) (cl:float b 1d0))))
199
200	#+cmu
201	(defmethod two-arg-/ ((a qd-real) (b ext:double-double-float))
202	(make-instance 'qd-real :value (div-qd-dd (qd-value a)
203	b)))
204
205	(defmethod two-arg-/ ((a cl:float) (b qd-real))
206	(make-instance 'qd-real :value (div-qd (make-qd-d (cl:float a 1d0))
207	(qd-value b))))
208
209	#+cmu
210	(defmethod two-arg-/ ((a ext:double-double-float) (b qd-real))
211	(make-instance 'qd-real :value (div-qd (make-qd-dd a 0w0)
212	(qd-value b))))
213
214	(defmethod two-arg-/ ((a number) (b number))
215	(cl:/ a b))
216
217	(defmethod unary-divide ((a number))
218	(cl:/ a))
219
220	(defmethod unary-divide ((a qd-real))
221	(make-instance 'qd-real :value (div-qd +qd-one+ (qd-value a))))
222
223	(defun / (number &rest more-numbers)
224	(if more-numbers
225	(do ((nlist more-numbers (cdr nlist))
226	(result number))
227	((atom nlist) result)
228	(declare (list nlist))
229	(setq result (two-arg-/ result (car nlist))))
230	(unary-divide number)))
231
232	(macrolet ((frob (name &optional (type 'real))
233	(let ((method-name (intern (concatenate 'string
234	(string '#:q)
235	(symbol-name name))))
236	(cl-name (intern (symbol-name name) :cl))
237	(qd-name (intern (concatenate 'string
238	(symbol-name name)
239	(string '#:-qd)))))
240	`(progn
241	(defmethod ,method-name ((x ,type))
242	(,cl-name x))
243	(defmethod ,method-name ((x qd-real))
244	(,qd-name (qd-value x)))
245	(declaim (inline ,name))
246	(defun ,name (x)
247	(,method-name x))))))
248	(frob zerop number)
249	(frob plusp)
250	(frob minusp))
251
252	(defun bignum-to-qd (bignum)
253	(make-instance 'qd-real
254	:value (rational-to-qd bignum)))
255
256	(defmethod qfloat ((x real) (num-type cl:float))
257	(cl:float x num-type))
258
259	(defmethod qfloat ((x cl:float) (num-type qd-real))
260	(make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
261
262	(defmethod qfloat ((x integer) (num-type qd-real))
263	(cond ((typep x 'fixnum)
264	(make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
265	(t
266	;; A bignum
267	(bignum-to-qd x))))
268
269	#+nil
270	(defmethod qfloat ((x ratio) (num-type qd-real))
271	;; This probably has some issues with roundoff
272	(two-arg-/ (qfloat (numerator x) num-type)
273	(qfloat (denominator x) num-type)))
274
275	(defmethod qfloat ((x ratio) (num-type qd-real))
276	(make-instance 'qd-real :value (rational-to-qd x)))
277
278	#+cmu
279	(defmethod qfloat ((x ext:double-double-float) (num-type qd-real))
280	(make-instance 'qd-real :value (make-qd-dd x 0w0)))
281
282	(defmethod qfloat ((x qd-real) (num-type cl:float))
283	(multiple-value-bind (q0 q1 q2 q3)
284	(qd-parts (qd-value x))
285	(cl:float (cl:+ q0 q1 q2 q3) num-type)))
286
287	#+cmu
288	(defmethod qfloat ((x qd-real) (num-type ext:double-double-float))
289	(multiple-value-bind (q0 q1 q2 q3)
290	(qd-parts (qd-value x))
291	(cl:+ (cl:float q0 1w0)
292	(cl:float q1 1w0)
293	(cl:float q2 1w0)
294	(cl:float q3 1w0))))
295
296	(defmethod qfloat ((x qd-real) (num-type qd-real))
297	x)
298
299	(declaim (inline float))
300	(defun float (x num-type)
301	(qfloat x num-type))
302
303	(defmethod qrealpart ((x number))
304	(cl:realpart x))
305	(defmethod qrealpart ((x qd-real))
306	x)
307	(defmethod qrealpart ((x qd-complex))
308	(make-instance 'qd-real :value (qd-real x)))
309	(defun realpart (x)
310	(qrealpart x))
311
312	(defmethod qimagpart ((x number))
313	(cl:imagpart x))
314	(defmethod qimagpart ((x qd-real))
315	(make-qd 0d0))
316	(defmethod qimagpart ((x qd-complex))
317	(make-instance 'qd-real :value (qd-imag x)))
318
319	(defun imagpart (x)
320	(qimagpart x))
321
322	(defmethod qconjugate ((a number))
323	(cl:conjugate a))
324
325	(defmethod qconjugate ((a qd-real))
326	(make-instance 'qd-real :value (qd-value a)))
327
328	(defmethod qconjugate ((a qd-complex))
329	(make-instance 'qd-complex
330	:real (qd-real a)
331	:imag (neg-qd (qd-imag a))))
332
333	(defun conjugate (z)
334	(qconjugate z))
335
336	(defmethod qscale-float ((f cl:float) (n integer))
337	(cl:scale-float f n))
338
339	(defmethod qscale-float ((f qd-real) (n integer))
340	(make-instance 'qd-real :value (scale-float-qd (qd-value f) n)))
341
342	(declaim (inline scale-float))
343	(defun scale-float (f n)
344	(qscale-float f n))
345
346	(macrolet
347	((frob (op)
348	(let ((method (intern (concatenate 'string
349	(string '#:two-arg-)
350	(symbol-name op))))
351	(cl-fun (find-symbol (symbol-name op) :cl))
352	(qd-fun (intern (concatenate 'string (string '#:qd-) (symbol-name op))
353	'#:octi)))
354	`(progn
355	(defmethod ,method ((a real) (b real))
356	(,cl-fun a b))
357	(defmethod ,method ((a qd-real) (b real))
358	(,qd-fun (qd-value a) (make-qd-d (cl:float b 1d0))))
359	(defmethod ,method ((a real) (b qd-real))
360	;; This is not really right if A is a rational. We're
361	;; supposed to compare them as rationals.
362	(,qd-fun (make-qd-d (cl:float a 1d0)) (qd-value b)))
363	(defmethod ,method ((a qd-real) (b qd-real))
364	(,qd-fun (qd-value a) (qd-value b)))
365	(defun ,op (number &rest more-numbers)
366	"Returns T if its arguments are in strictly increasing order, NIL otherwise."
367	(declare (optimize (safety 2))
368	(dynamic-extent more-numbers))
369	(do* ((n number (car nlist))
370	(nlist more-numbers (cdr nlist)))
371	((atom nlist) t)
372	(declare (list nlist))
373	(if (not (,method n (car nlist))) (return nil))))))))
374	(frob <)
375	(frob >)
376	(frob <=)
377	(frob >=))
378
379	;; Handle the special functions for a real argument. Complex args are
380	;; handled elsewhere.
381	(macrolet
382	((frob (name)
383	(let ((method-name
384	(intern (concatenate 'string (string '#:q)
385	(symbol-name name))))
386	(cl-name (intern (symbol-name name) :cl))
387	(qd-name (intern (concatenate 'string (symbol-name name)
388	(string '#:-qd)))))
389	`(progn
390	(defmethod ,name ((x number))
391	(,cl-name x))
392	(defmethod ,name ((x qd-real))
393	(make-instance 'qd-real :value (,qd-name (qd-value x))))))))
394	(frob abs)
395	(frob exp)
396	(frob sin)
397	(frob cos)
398	(frob tan)
399	;;(frob asin)
400	;;(frob acos)
401	(frob sinh)
402	(frob cosh)
403	(frob tanh)
404	(frob asinh)
405	;;(frob acosh)
406	;;(frob atanh)
407	)
408
409	(defmethod sqrt ((x number))
410	(cl:sqrt x))
411
412	(defmethod sqrt ((x qd-real))
413	(if (minusp x)
414	(make-instance 'qd-complex
415	:real +qd-zero+
416	:imag (sqrt-qd (neg-qd (qd-value x))))
417	(make-instance 'qd-real :value (sqrt-qd (qd-value x)))))
418
419	(defun scalb (x n)
420	"Compute 2^N * X without compute 2^N first (use properties of the
421	underlying floating-point format"
422	(declare (type qd-real x))
423	(scale-float x n))
424
425	(declaim (inline qd-cssqs))
426	(defun qd-cssqs (z)
427	(multiple-value-bind (rho k)
428	(octi::hypot-aux-qd (qd-value (realpart z))
429	(qd-value (imagpart z)))
430	(values (make-instance 'qd-real :value rho)
431	k)))
432
433	#+nil
434	(defmethod qabs ((z qd-complex))
435	;; sqrt(x^2+y^2)
436	;; If \|x\| > \|y\| then sqrt(x^2+y^2) = \|x\|*sqrt(1+(y/x)^2)
437	(multiple-value-bind (abs^2 rho)
438	(hypot-qd (qd-value (realpart z))
439	(qd-value (imagpart z)))
440	(scale-float (make-instance 'qd-real :value (sqrt abs^2))
441	rho)))
442
443	(defmethod abs ((z qd-complex))
444	;; sqrt(x^2+y^2)
445	;; If \|x\| > \|y\| then sqrt(x^2+y^2) = \|x\|*sqrt(1+(y/x)^2)
446	(make-instance 'qd-real
447	:value (hypot-qd (qd-value (realpart z))
448	(qd-value (imagpart z)))))
449
450	(defmethod log ((a number) &optional b)
451	(if b
452	(cl:log a b)
453	(cl:log a)))
454
455	(defmethod log ((a qd-real) &optional b)
456	(if b
457	(/ (log a) (log b))
458	(if (minusp (float-sign a))
459	(make-instance 'qd-complex
460	:real (log-qd (abs-qd (qd-value a)))
461	:imag +qd-pi+)
462	(make-instance 'qd-real :value (log-qd (qd-value a))))))
463
464	(defmethod log1p ((a qd-real))
465	(make-instance 'qd-real :value (log1p-qd (qd-value a))))
466
467	(defmethod atan ((y real) &optional x)
468	(cond (x
469	(cond ((typep x 'qd-real)
470	(make-instance 'qd-real
471	:value (atan2-qd (qd-value y) (qd-value x))))
472	(t
473	(cl:atan y x))))
474	(t
475	(cl:atan y))))
476
477	(defmethod atan ((y qd-real) &optional x)
478	(make-instance 'qd-real
479	:value
480	(if x
481	(atan2-qd (qd-value y) (qd-value x))
482	(atan-qd (qd-value y)))))
483
484	(defmethod qexpt ((x number) (y number))
485	(cl:expt x y))
486
487	(defmethod qexpt ((x qd-real) (y real))
488	(exp (* y (log x))))
489
490	(defmethod qexpt ((x real) (y qd-real))
491	(exp (* y (log x))))
492
493	(defmethod qexpt ((x qd-real) (y cl:complex))
494	(exp (* (make-instance 'qd-complex
495	:real (qd-value (realpart y))
496	:imag (qd-value (imagpart y)))
497	(log x))))
498
499	(defmethod qexpt ((x cl:complex) (y qd-real))
500	(exp (* y
501	(log (make-instance 'qd-complex
502	:real (qd-value (realpart x))
503	:imag (qd-value (imagpart x)))))))
504
505	(defmethod qexpt ((x qd-real) (y qd-real))
506	;; x^y = exp(y*log(x))
507	(exp (* y (log x))))
508
509	(defmethod qexpt ((x qd-real) (y integer))
510	(make-instance 'qd-real
511	:value (npow (qd-value x) y)))
512
513	(declaim (inline expt))
514	(defun expt (x y)
515	(qexpt x y))
516
517
518
519	(defmethod two-arg-= ((a number) (b number))
520	(cl:= a b))
521
522	(defmethod two-arg-= ((a qd-real) (b number))
523	(if (cl:realp b)
524	(qd-= (qd-value a) (make-qd-d (cl:float b 1d0)))
525	nil))
526
527	(defmethod two-arg-= ((a number) (b qd-real))
528	(if (cl:realp a)
529	(qd-= (make-qd-d (cl:float a 1d0)) (qd-value b))
530	nil))
531
532	(defmethod two-arg-= ((a qd-complex) b)
533	(and (two-arg-= (realpart a) (realpart b))
534	(two-arg-= (imagpart a) (imagpart b))))
535
536	(defmethod two-arg-= (a (b qd-complex))
537	(and (two-arg-= (realpart a) (realpart b))
538	(two-arg-= (imagpart a) (imagpart b))))
539
540
541	(defmethod two-arg-= ((a qd-real) (b qd-real))
542	(qd-= (qd-value a) (qd-value b)))
543
544	(defun = (number &rest more-numbers)
545	"Returns T if all of its arguments are numerically equal, NIL otherwise."
546	(declare (optimize (safety 2))
547	(dynamic-extent more-numbers))
548	(do ((nlist more-numbers (cdr nlist)))
549	((atom nlist) t)
550	(declare (list nlist))
551	(if (not (two-arg-= (car nlist) number))
552	(return nil))))
553
554	(defun /= (number &rest more-numbers)
555	"Returns T if no two of its arguments are numerically equal, NIL otherwise."
556	(declare (optimize (safety 2))
557	(dynamic-extent more-numbers))
558	(do* ((head number (car nlist))
559	(nlist more-numbers (cdr nlist)))
560	((atom nlist) t)
561	(declare (list nlist))
562	(unless (do* ((nl nlist (cdr nl)))
563	((atom nl) t)
564	(declare (list nl))
565	(if (two-arg-= head (car nl))
566	(return nil)))
567	(return nil))))
568
569	(defmethod qcomplex ((x real) &optional y)
570	(cl:complex x (if y y 0)))
571
572	(defmethod qcomplex ((x qd-real) &optional y)
573	(make-instance 'qd-complex
574	:real (qd-value x)
575	:imag (if y (qd-value y) +qd-zero+)))
576
577	(defun complex (x &optional (y 0))
578	(qcomplex x y))
579
580	(defmethod qinteger-decode-float ((f cl:float))
581	(cl:integer-decode-float f))
582
583	(defmethod qinteger-decode-float ((f qd-real))
584	(integer-decode-qd (qd-value f)))
585
586	(declaim (inline integer-decode-float))
587	(defun integer-decode-float (f)
588	(qinteger-decode-float f))
589
590	(defmethod qdecode-float ((f cl:float))
591	(cl:decode-float f))
592
593	(defmethod qdecode-float ((f qd-real))
594	(multiple-value-bind (frac exp s)
595	(decode-float-qd (qd-value f))
596	(values (make-instance 'qd-real :value frac)
597	exp
598	(make-instance 'qd-real :value s))))
599
600	(declaim (inline decode-float))
601	(defun decode-float (f)
602	(qdecode-float f))
603
604	(defmethod qfloor ((x real) &optional y)
605	(if y
606	(cl:floor x y)
607	(cl:floor x)))
608
609	(defmethod qfloor ((x qd-real) &optional y)
610	(if (and y (/= y 1))
611	(let ((f (qfloor (/ x y))))
612	(values f
613	(- x (* f y))))
614	(let ((f (ffloor-qd (qd-value x))))
615	(multiple-value-bind (int exp sign)
616	(integer-decode-qd f)
617	(values (ash (* sign int) exp)
618	(make-instance 'qd-real
619	:value (qd-value
620	(- x (make-instance 'qd-real
621	:value f)))))))))
622
623	(defun floor (x &optional y)
624	(qfloor x y))
625
626	(defmethod qffloor ((x real) &optional y)
627	(if y
628	(cl:ffloor x y)
629	(cl:ffloor x)))
630
631	(defmethod qffloor ((x qd-real) &optional y)
632	(if (and y (/= y 1))
633	(let ((f (qffloor (/ x y))))
634	(values f
635	(- x (* f y))))
636	(let ((f (make-instance 'qd-real :value (ffloor-qd (qd-value x)))))
637	(values f
638	(- x f)))))
639
640	(defun ffloor (x &optional y)
641	(qffloor x y))
642
643	(defun ceiling (x &optional y)
644	(multiple-value-bind (f rem)
645	(floor x y)
646	(if (zerop rem)
647	(values (+ f 1)
648	rem)
649	(values (+ f 1)
650	(- rem 1)))))
651
652	(defun fceiling (x &optional y)
653	(multiple-value-bind (f rem)
654	(ffloor x y)
655	(if (zerop rem)
656	(values (+ f 1)
657	rem)
658	(values (+ f 1)
659	(- rem 1)))))
660
661	(defun truncate (x &optional (y 1))
662	(if (minusp x)
663	(ceiling x y)
664	(floor x y)))
665
666	(defun ftruncate (x &optional (y 1))
667	(if (minusp x)
668	(fceiling x y)
669	(ffloor x y)))
670
671	(defmethod %unary-round ((x real))
672	(cl::round x))
673
674	(defmethod %unary-round ((number qd-real))
675	(multiple-value-bind (bits exp)
676	(integer-decode-float number)
677	(let* ((shifted (ash bits exp))
678	(rounded (if (and (minusp exp)
679	(oddp shifted)
680	(not (zerop (logand bits
681	(ash 1 (- -1 exp))))))
682	(1+ shifted)
683	shifted)))
684	(if (minusp number)
685	(- rounded)
686	rounded))))
687
688	(defun round (number &optional (divisor 1))
689	(if (eql divisor 1)
690	(let ((r (%unary-round number)))
691	(values r
692	(- number r)))
693	(multiple-value-bind (tru rem)
694	(truncate number divisor)
695	(if (zerop rem)
696	(values tru rem)
697	(let ((thresh (/ (abs divisor) 2)))
698	(cond ((or (> rem thresh)
699	(and (= rem thresh) (oddp tru)))
700	(if (minusp divisor)
701	(values (- tru 1) (+ rem divisor))
702	(values (+ tru 1) (- rem divisor))))
703	((let ((-thresh (- thresh)))
704	(or (< rem -thresh)
705	(and (= rem -thresh) (oddp tru))))
706	(if (minusp divisor)
707	(values (+ tru 1) (- rem divisor))
708	(values (- tru 1) (+ rem divisor))))
709	(t (values tru rem))))))))
710
711	(defun fround (number &optional (divisor 1))
712	"Same as ROUND, but returns first value as a float."
713	(multiple-value-bind (res rem)
714	(round number divisor)
715	(values (float res (if (floatp rem) rem 1.0)) rem)))
716
717	(defmethod qfloat-sign ((a real) &optional (f (float 1 a)))
718	(cl:float-sign a f))
719
720
721	(defmethod qfloat-sign ((a qd-real) &optional f)
722	(if f
723	(make-instance 'qd-real
724	:value (mul-qd-d (abs-qd (qd-value f))
725	(cl:float-sign (qd-0 (qd-value a)))))
726	(make-instance 'qd-real :value (make-qd-d (cl:float-sign (qd-0 (qd-value a)))))))
727
728	(declaim (inline float-sign))
729	(defun float-sign (n &optional (float2 nil float2p))
730	(if float2p
731	(qfloat-sign n float2)
732	(qfloat-sign n)))
733
734	(defun max (number &rest more-numbers)
735	"Returns the greatest of its arguments."
736	(declare (optimize (safety 2)) (type (or real qd-real) number)
737	(dynamic-extent more-numbers))
738	(dolist (real more-numbers)
739	(when (> real number)
740	(setq number real)))
741	number)
742
743	(defun min (number &rest more-numbers)
744	"Returns the least of its arguments."
745	(declare (optimize (safety 2)) (type (or real qd-real) number)
746	(dynamic-extent more-numbers))
747	(do ((nlist more-numbers (cdr nlist))
748	(result (the (or real qd-real) number)))
749	((null nlist) (return result))
750	(declare (list nlist))
751	(if (< (car nlist) result)
752	(setq result (car nlist)))))
753
754	(defmethod asin ((x number))
755	(cl:asin x))
756
757	(defmethod asin ((x qd-real))
758	(if (<= -1 x 1)
759	(make-instance 'qd-real :value (asin-qd (qd-value x)))
760	(qd-complex-asin x)))
761
762	(defmethod acos ((x number))
763	(cl:acos x))
764
765	(defmethod acos ((x qd-real))
766	(cond ((> (abs x) 1)
767	(qd-complex-acos x))
768	(t
769	(make-instance 'qd-real :value (acos-qd (qd-value x))))))
770
771	(defmethod acosh ((x number))
772	(cl:acosh x))
773
774	(defmethod acosh ((x qd-real))
775	(if (< x 1)
776	(qd-complex-acosh x)
777	(make-instance 'qd-real :value (acosh-qd (qd-value x)))))
778
779	(defmethod atanh ((x number))
780	(cl:atanh x))
781
782	(defmethod atanh ((x qd-real))
783	(if (> (abs x) 1)
784	(qd-complex-atanh x)
785	(make-instance 'qd-real :value (atanh-qd (qd-value x)))))
786
787	(defmethod cis ((x real))
788	(cl:cis x))
789
790	(defmethod cis ((x qd-real))
791	(multiple-value-bind (s c)
792	(sincos-qd (qd-value x))
793	(make-instance 'qd-complex
794	:real c
795	:imag s)))
796
797	(defmethod phase ((x number))
798	(cl:phase x))
799
800	(defmethod phase ((x qd-real))
801	(if (minusp x)
802	(- +pi+)
803	(make-instance 'qd-real :value (make-qd-d 0d0))))
804
805	(defun signum (number)
806	"If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
807	(if (zerop number)
808	number
809	(if (rationalp number)
810	(if (plusp number) 1 -1)
811	(/ number (abs number)))))
812
813	(defmethod random ((x cl:real) &optional (state random-state))
814	(cl:random x state))
815
816	(defmethod random ((x qd-real) &optional (state random-state))
817	(* x (make-instance 'qd-real
818	:value (octi:random-qd state))))
819
820	(defmethod float-digits ((x cl:real))
821	(cl:float-digits x))
822
823	(defmethod float-digits ((x qd-real))
824	(* 4 (float-digits 1d0)))
825
826	(defmethod rational ((x real))
827	(cl:rational x))
828
829	(defmethod rational ((x qd-real))
830	(with-qd-parts (x0 x1 x2 x3)
831	(qd-value x)
832	(+ (cl:rational x0)
833	(cl:rational x1)
834	(cl:rational x2)
835	(cl:rational x3))))
836
837	(defmethod rationalize ((x cl:real))
838	(cl:rationalize x))
839
840	;;; The algorithm here is the method described in CLISP. Bruno Haible has
841	;;; graciously given permission to use this algorithm. He says, "You can use
842	;;; it, if you present the following explanation of the algorithm."
843	;;;
844	;;; Algorithm (recursively presented):
845	;;; If x is a rational number, return x.
846	;;; If x = 0.0, return 0.
847	;;; If x < 0.0, return (- (rationalize (- x))).
848	;;; If x > 0.0:
849	;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
850	;;; exponent, sign).
851	;;; If m = 0 or e >= 0: return x = m*2^e.
852	;;; Search a rational number between a = (m-1/2)2^e and b = (m+1/2)2^e
853	;;; with smallest possible numerator and denominator.
854	;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
855	;;; But in this case the result will be x itself anyway, regardless of
856	;;; the choice of a. Therefore we can simply ignore this case.
857	;;; Note 2: At first, we need to consider the closed interval [a,b].
858	;;; but since a and b have the denominator 2^(\|e\|+1) whereas x itself
859	;;; has a denominator <= 2^\|e\|, we can restrict the seach to the open
860	;;; interval (a,b).
861	;;; So, for given a and b (0 < a < b) we are searching a rational number
862	;;; y with a <= y <= b.
863	;;; Recursive algorithm fraction_between(a,b):
864	;;; c := (ceiling a)
865	;;; if c < b
866	;;; then return c ; because a <= c < b, c integer
867	;;; else
868	;;; ; a is not integer (otherwise we would have had c = a < b)
869	;;; k := c-1 ; k = floor(a), k < a < b <= k+1
870	;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
871	;;; ; note 1 <= 1/(b-k) < 1/(a-k)
872	;;;
873	;;; You can see that we are actually computing a continued fraction expansion.
874	;;;
875	;;; Algorithm (iterative):
876	;;; If x is rational, return x.
877	;;; Call (integer-decode-float x). It returns a m,e,s (mantissa,
878	;;; exponent, sign).
879	;;; If m = 0 or e >= 0, return m2^es. (This includes the case x = 0.0.)
880	;;; Create rational numbers a := (2m-1)2^(e-1) and b := (2m+1)2^(e-1)
881	;;; (positive and already in lowest terms because the denominator is a
882	;;; power of two and the numerator is odd).
883	;;; Start a continued fraction expansion
884	;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
885	;;; Loop
886	;;; c := (ceiling a)
887	;;; if c >= b
888	;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
889	;;; goto Loop
890	;;; finally partial_quotient(c).
891	;;; Here partial_quotient(c) denotes the iteration
892	;;; i := i+1, p[i] := cp[i-1]+p[i-2], q[i] := cq[i-1]+q[i-2].
893	;;; At the end, return s * (p[i]/q[i]).
894	;;; This rational number is already in lowest terms because
895	;;; p[i]q[i-1]-p[i-1]q[i] = (-1)^i.
896	;;;
897	(defmethod rationalize ((x qd-real))
898	;; This is a fairly straigtforward implementation of the iterative
899	;; algorithm above.
900	(multiple-value-bind (frac expo sign)
901	(integer-decode-float x)
902	(cond ((or (zerop frac) (>= expo 0))
903	(if (minusp sign)
904	(- (ash frac expo))
905	(ash frac expo)))
906	(t
907	;; expo < 0 and (2m-1) and (2m+1) are coprime to 2^(1-e),
908	;; so build the fraction up immediately, without having to do
909	;; a gcd.
910	(let ((a (/ (- (* 2 frac) 1) (ash 1 (- 1 expo))))
911	(b (/ (+ (* 2 frac) 1) (ash 1 (- 1 expo))))
912	(p0 0)
913	(q0 1)
914	(p1 1)
915	(q1 0))
916	(do ((c (ceiling a) (ceiling a)))
917	((< c b)
918	(let ((top (+ (* c p1) p0))
919	(bot (+ (* c q1) q0)))
920	(/ (if (minusp sign)
921	(- top)
922	top)
923	bot)))
924	(let* ((k (- c 1))
925	(p2 (+ (* k p1) p0))
926	(q2 (+ (* k q1) q0)))
927	(psetf a (/ (- b k))
928	b (/ (- a k)))
929	(setf p0 p1
930	q0 q1
931	p1 p2
932	q1 q2))))))))
933
934	(define-compiler-macro + (&whole form &rest args)
935	(declare (ignore form))
936	(if (null args)
937	0
938	(do ((args (cdr args) (cdr args))
939	(res (car args)
940	`(two-arg-+ ,res ,(car args))))
941	((null args) res))))
942
943	(define-compiler-macro - (&whole form number &rest more-numbers)
944	(declare (ignore form))
945	(if more-numbers
946	(do ((nlist more-numbers (cdr nlist))
947	(result number))
948	((atom nlist) result)
949	(declare (list nlist))
950	(setq result `(two-arg-- ,result ,(car nlist))))
951	`(unary-minus ,number)))
952
953	(define-compiler-macro * (&whole form &rest args)
954	(declare (ignore form))
955	(if (null args)
956	1
957	(do ((args (cdr args) (cdr args))
958	(res (car args)
959	`(two-arg-* ,res ,(car args))))
960	((null args) res))))
961
962	(define-compiler-macro / (number &rest more-numbers)
963	(if more-numbers
964	(do ((nlist more-numbers (cdr nlist))
965	(result number))
966	((atom nlist) result)
967	(declare (list nlist))
968	(setq result `(two-arg-/ ,result ,(car nlist))))
969	`(unary-divide ,number)))
970
971	;; Compiler macros to convert <, >, <=, and >= into multiple calls of
972	;; the corresponding two-arg-<foo> function.
973	(macrolet
974	((frob (op)
975	(let ((method (intern (concatenate 'string
976	(string '#:two-arg-)
977	(symbol-name op)))))
978	`(define-compiler-macro ,op (number &rest more-numbers)
979	(do* ((n number (car nlist))
980	(nlist more-numbers (cdr nlist))
981	(res nil))
982	((atom nlist)
983	`(and ,@(nreverse res)))
984	(push `(,',method ,n ,(car nlist)) res))))))
985	(frob <)
986	(frob >)
987	(frob <=)
988	(frob >=))
989
990	(define-compiler-macro /= (&whole form number &rest more-numbers)
991	;; Convert (/= x y) to (not (two-arg-= x y)). Should we try to
992	;; handle a few more cases?
993	(if (cdr more-numbers)
994	form
995	`(not (two-arg-= ,number ,(car more-numbers)))))
996
997
998	;; Define compiler macro the convert two-arg-foo into the appropriate
999	;; CL function or QD-REAL function so we don't have to do CLOS
1000	;; dispatch.
1001	#+(or)
1002	(macrolet
1003	((frob (name cl-op qd-op)
1004	`(define-compiler-macro ,name (&whole form x y &environment env)
1005	(flet ((arg-type (arg)
1006	(multiple-value-bind (def-type localp decl)
1007	(ext:variable-information arg env)
1008	(declare (ignore localp))
1009	(when def-type
1010	(cdr (assoc 'type decl))))))
1011	(let ((x-type (arg-type x))
1012	(y-type (arg-type y)))
1013	(cond ((and (subtypep x-type 'cl:number)
1014	(subtypep y-type 'cl:number))
1015	`(,',cl-op ,x ,y))
1016	((and (subtypep x-type 'qd-real)
1017	(subtypep y-type 'qd-real))
1018	`(make-instance 'qd-real :value (,',qd-op (qd-value ,x)
1019	(qd-value ,y))))
1020	(t
1021	;; Don't know how to handle this, so give up.
1022	form)))))))
1023	(frob two-arg-+ cl:+ add-qd)
1024	(frob two-arg-- cl:- sub-qd)
1025	(frob two-arg-* cl:* mul-qd)
1026	(frob two-arg-/ cl:/ div-qd))
1027
1028	#+(or)
1029	(macrolet
1030	((frob (name cl-op qd-op cl-qd-op qd-cl-op)
1031	`(define-compiler-macro ,name (&whole form x y &environment env)
1032	(flet ((arg-type (arg)
1033	(multiple-value-bind (def-type localp decl)
1034	(ext:variable-information arg env)
1035	(declare (ignore localp))
1036	(when def-type
1037	(cdr (assoc 'type decl))))))
1038	(let ((x-type (arg-type x))
1039	(y-type (arg-type y)))
1040	(cond ((subtypep x-type 'cl:float)
1041	(cond ((subtypep y-type 'cl:number)
1042	`(,',cl-op ,x ,y))
1043	((subtypep y-type 'qd-real)
1044	(if ,cl-qd-op
1045	`(make-instance 'qd-real :value (,',cl-qd-op (cl:float ,x 1d0)
1046	(qd-value ,y)))
1047	form))
1048	(t form)))
1049	((subtypep x-type 'qd-real)
1050	(cond ((subtypep y-type 'cl:float)
1051	(if ,qd-cl-op
1052	`(make-instance 'qd-real :value (,',qd-cl-op (qd-value ,x)
1053	(float ,y 1d0)))
1054	form))
1055	((subtypep y-type 'qd-real)
1056	`(make-instance 'qd-real :value (,',qd-op (qd-value ,x)
1057	(qd-value ,y))))
1058	(t form)))
1059	(t
1060	;; Don't know how to handle this, so give up.
1061	form)))))))
1062	(frob two-arg-+ cl:+ add-qd add-d-qd add-qd-d)
1063	(frob two-arg-- cl:- sub-qd sub-d-qd sub-qd-d)
1064	(frob two-arg-* cl:* mul-qd mul-d-qd mul-qd-d)
1065	(frob two-arg-/ cl:/ div-qd nil nil))
1066
1067
1068	(defun read-qd-real-or-complex (stream)
1069	(let ((c (peek-char t stream)))
1070	(cond ((char= c #\()
1071	;; Read a QD complex
1072	(read-char stream) ; Skip the paren
1073	(let ((real (read stream t nil t))
1074	(imag (read stream t nil t)))
1075	(unless (char= (peek-char t stream) #\))
1076	(error "Illegal QD-COMPLEX number format"))
1077	;; Read closing paren
1078	(read-char stream)
1079	(make-instance 'qd-complex
1080	:real (qd-value (float real +qd-real-one+))
1081	:imag (qd-value (float imag +qd-real-one+)))))
1082	(t
1083	(make-instance 'qd-real :value (read-qd stream))))))
1084
1085	(defun qd-class-reader (stream subchar arg)
1086	(declare (ignore subchar))
1087	(when arg
1088	(warn "Numeric argument ignored in #~DQ" arg))
1089	(read-qd-real-or-complex stream))
1090
1091	;; Yow! We redefine the #q reader that is in qd-io.lisp to read in
1092	;; and make a real qd-real float, instead of the hackish
1093	;; %qd-real.
1094	(set-dispatch-macro-character #\# #\Q #'qd-class-reader)
1095
1096
1097	(defmethod epsilon ((m cl:float))
1098	(etypecase m
1099	(single-float single-float-epsilon)
1100	(double-float double-float-epsilon)))
1101
1102	(defmethod epsilon ((m cl:complex))
1103	(epsilon (realpart m)))
1104
1105	(defmethod epsilon ((m qd-real))
1106	;; What is the epsilon value for a quad-double? This is complicated
1107	;; by the fact that things like (+ #q1 #q1q-100) is representable as
1108	;; a quad-double. For most purposes we want epsilon to be close to
1109	;; the 212 bits of precision (4*53 bits) that we normally have with
1110	;; a quad-double.
1111	(scale-float (make-qd-d 1d0) -212))
1112
1113	(defmethod epsilon ((m qd-complex))
1114	(epsilon (realpart m)))

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format