root/qd-methods.lisp @ 1e343ca0be49b9e56938a8063b4d2e679e338965

Revision 1e343ca0be49b9e56938a8063b4d2e679e338965, 32.0 KB (checked in by Raymond Toy <toy.raymond@…>, 3 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
421underlying 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 m*2^e*s. (This includes the case x = 0.0.)
880;;;   Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+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] := c*p[i-1]+p[i-2], q[i] := c*q[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 (2*m-1) and (2*m+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.