Context Navigation

root/qd-methods.lisp @ bd6d814332fdf6a831873bb33c9e4753f29d414f

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

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format