source: qd-methods.lisp @ bd6d81

Last change on this file since bd6d81 was bd6d81, checked in by Raymond Toy <toy.raymond@…>, 2 years ago

Oops. The second arg to FLOAT is optional!

  • Property mode set to 100644
File size: 32.2 KB
Line 
1;;;; -*- Mode: lisp -*-
2;;;;
3;;;; Copyright (c) 2007, 2008, 2011 Raymond Toy
4;;;;
5;;;; Permission is hereby granted, free of charge, to any person
6;;;; obtaining a copy of this software and associated documentation
7;;;; files (the "Software"), to deal in the Software without
8;;;; restriction, including without limitation the rights to use,
9;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10;;;; copies of the Software, and to permit persons to whom the
11;;;; Software is furnished to do so, subject to the following
12;;;; conditions:
13;;;;
14;;;; The above copyright notice and this permission notice shall be
15;;;; included in all copies or substantial portions of the Software.
16;;;;
17;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24;;;; OTHER DEALINGS IN THE SOFTWARE.
25
26(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
423underlying 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 m*2^e*s. (This includes the case x = 0.0.)
893;;;   Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+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] := c*p[i-1]+p[i-2], q[i] := c*q[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 (2*m-1) and (2*m+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
1083the smallest number x such that 1+x /= x.  The argument can be
1084complex"))
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
1108the 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 repository browser.