external/oct/qd-methods.lisp

   1 ;;;; -*- Mode: lisp -*-
   2 ;;;;
   3 ;;;; Copyright (c) 2007 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 #:qd)
  27
  28 (defconstant +pi+
  29   (make-instance 'qd-real :value qdi:+qd-pi+))
  30
  31 #+cmu
  32 (defconstant +quad-double-float-positive-infinity+
  33   (make-instance 'qd-real :value (make-qd-d ext:double-float-positive-infinity))
  34   "Positive infinity for qd-real")
  35
  36 #+cmu
  37 (defconstant +quad-double-float-negative-infinity+
  38   (make-instance 'qd-real :value (make-qd-d ext:double-float-negative-infinity))
  39   "Negative infinity for qd-real")
  40
  41 (defconstant +most-positive-quad-double-float+
  42   (make-instance 'qd-real
  43                  :value (qdi::%make-qd-d most-positive-double-float
  44                                          (cl:scale-float most-positive-double-float (cl:* 1 -53))
  45                                          (cl:scale-float most-positive-double-float (cl:* 2 -53))
  46                                          (cl:scale-float most-positive-double-float (cl:* 3 -53)))))
  47
  48 (defconstant +least-positive-quad-double-float+
  49   (make-instance 'qd-real
  50                  :value (make-qd-d least-positive-double-float)))
  51
  52 ;; Not sure this is 100% correct, but I think if the first component
  53 ;; is any smaller than this, the last component would no longer be a
  54 ;; normalized double-float.
  55 (defconstant +least-positive-normalized-quad-double-float+
  56   (make-instance 'qd-real
  57                  :value (make-qd-d (cl:scale-float least-positive-normalized-double-float (cl:* 3 53)))))
  58
  59 (defconstant +qd-real-one+
  60   (make-instance 'qd-real :value (make-qd-d 1d0)))
  61
  62 (defmethod add1 ((a number))
  63   (cl::1+ a))
  64
  65 (defmethod add1 ((a qd-real))
  66   (make-instance 'qd-real :value (add-qd-d (qd-value a) 1d0)))
  67
  68 (defmethod sub1 ((a number))
  69   (cl::1- a))
  70
  71 (defmethod sub1 ((a qd-real))
  72   (make-instance 'qd-real :value (sub-qd-d (qd-value a) 1d0)))
  73
  74 (declaim (inline 1+ 1-))
  75
  76 (defun 1+ (x)
  77   (add1 x))
  78
  79 (defun 1- (x)
  80   (sub1 x))
  81
  82 (defmethod two-arg-+ ((a qd-real) (b qd-real))
  83   (make-instance 'qd-real :value (add-qd (qd-value a) (qd-value b))))
  84
  85 (defmethod two-arg-+ ((a qd-real) (b cl:float))
  86   (make-instance 'qd-real :value (add-qd-d (qd-value a) (cl:float b 1d0))))
  87
  88 #+cmu
  89 (defmethod two-arg-+ ((a qd-real) (b ext:double-double-float))
  90   (make-instance 'qd-real :value (add-qd-dd (qd-value a) b)))
  91
  92 (defmethod two-arg-+ ((a real) (b qd-real))
  93   (+ b a))
  94
  95 (defmethod two-arg-+ ((a number) (b number))
  96   (cl:+ a b))
  97
  98 (defun + (&rest args)
  99   (if (null args)
 100       0
 101       (do ((args (cdr args) (cdr args))
 102            (res (car args)
 103                 (two-arg-+ res (car args))))
 104           ((null args) res))))
 105
 106 (defmethod two-arg-- ((a qd-real) (b qd-real))
 107   (make-instance 'qd-real :value (sub-qd (qd-value a) (qd-value b))))
 108
 109 (defmethod two-arg-- ((a qd-real) (b cl:float))
 110   (make-instance 'qd-real :value (sub-qd-d (qd-value a) (cl:float b 1d0))))
 111
 112 #+cmu
 113 (defmethod two-arg-- ((a qd-real) (b ext:double-double-float))
 114   (make-instance 'qd-real :value (sub-qd-dd (qd-value a) b)))
 115
 116 (defmethod two-arg-- ((a cl:float) (b qd-real))
 117   (make-instance 'qd-real :value (sub-d-qd (cl:float a 1d0) (qd-value b))))
 118
 119 (defmethod two-arg-- ((a number) (b number))
 120   (cl:- a b))
 121
 122 (defmethod unary-minus ((a number))
 123   (cl:- a))
 124
 125 (defmethod unary-minus ((a qd-real))
 126   (make-instance 'qd-real :value (neg-qd (qd-value a))))
 127
 128 (defun - (number &rest more-numbers)
 129   (if more-numbers
 130       (do ((nlist more-numbers (cdr nlist))
 131            (result number))
 132           ((atom nlist) result)
 133          (declare (list nlist))
 134          (setq result (two-arg-- result (car nlist))))
 135       (unary-minus number)))
 136
 137
 138 (defmethod two-arg-* ((a qd-real) (b qd-real))
 139   (make-instance 'qd-real :value (mul-qd (qd-value a) (qd-value b))))
 140
 141 (defmethod two-arg-* ((a qd-real) (b cl:float))
 142   (make-instance 'qd-real :value (mul-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   ;; We'd normally want to use mul-qd-dd, but mul-qd-dd is broken.
 147   (make-instance 'qd-real :value (mul-qd (qd-value a)
 148                                          (make-qd-dd b 0w0))))
 149
 150 (defmethod two-arg-* ((a real) (b qd-real))
 151   (* b a))
 152
 153 (defmethod two-arg-* ((a number) (b number))
 154   (cl:* a b))
 155
 156 (defun * (&rest args)
 157   (if (null args)
 158       1
 159       (do ((args (cdr args) (cdr args))
 160            (res (car args)
 161                 (two-arg-* res (car args))))
 162           ((null args) res))))
 163
 164 (defmethod two-arg-/ ((a qd-real) (b qd-real))
 165   (make-instance 'qd-real :value (div-qd (qd-value a) (qd-value b))))
 166
 167 (defmethod two-arg-/ ((a qd-real) (b cl:float))
 168   (make-instance 'qd-real :value (div-qd-d (qd-value a) (cl:float b 1d0))))
 169
 170 #+cmu
 171 (defmethod two-arg-/ ((a qd-real) (b ext:double-double-float))
 172   (make-instance 'qd-real :value (div-qd-dd (qd-value a)
 173                                             b)))
 174
 175 (defmethod two-arg-/ ((a cl:float) (b qd-real))
 176   (make-instance 'qd-real :value (div-qd (make-qd-d (cl:float a 1d0))
 177                                          (qd-value b))))
 178
 179 #+cmu
 180 (defmethod two-arg-/ ((a ext:double-double-float) (b qd-real))
 181   (make-instance 'qd-real :value (div-qd (make-qd-dd a 0w0)
 182                                          (qd-value b))))
 183
 184 (defmethod two-arg-/ ((a number) (b number))
 185   (cl:/ a b))
 186
 187 (defmethod unary-divide ((a number))
 188   (cl:/ a))
 189
 190 (defmethod unary-divide ((a qd-real))
 191   (make-instance 'qd-real :value (div-qd +qd-one+ (qd-value a))))
 192
 193 (defun / (number &rest more-numbers)
 194   (if more-numbers
 195       (do ((nlist more-numbers (cdr nlist))
 196            (result number))
 197           ((atom nlist) result)
 198          (declare (list nlist))
 199          (setq result (two-arg-/ result (car nlist))))
 200       (unary-divide number)))
 201
 202 (macrolet ((frob (name &optional (type 'real))
 203              (let ((method-name (intern (concatenate 'string "Q" (symbol-name name))))
 204                    (cl-name (intern (symbol-name name) :cl))
 205                    (qd-name (intern (concatenate 'string (symbol-name name) "-QD"))))
 206                `(progn
 207                   (defmethod ,method-name ((x ,type))
 208                     (,cl-name x))
 209                   (defmethod ,method-name ((x qd-real))
 210                     (,qd-name (qd-value x)))
 211                   (declaim (inline ,name))
 212                   (defun ,name (x)
 213                     (,method-name x))))))
 214   (frob zerop number)
 215   (frob plusp)
 216   (frob minusp))
 217
 218 (defun bignum-to-qd (bignum)
 219   (make-instance 'qd-real
 220                  :value (qdi::make-float (if (minusp bignum) -1 1)
 221                                          (abs bignum)
 222                                          0
 223                                          0
 224                                          0)))
 225
 226 (defmethod qfloat ((x real) (num-type cl:float))
 227   (cl:float x num-type))
 228
 229 (defmethod qfloat ((x cl:float) (num-type qd-real))
 230   (make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
 231
 232 (defmethod qfloat ((x integer) (num-type qd-real))
 233   (cond ((typep x 'fixnum)
 234          (make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
 235         (t
 236          ;; A bignum
 237          (bignum-to-qd x))))
 238
 239 #+nil
 240 (defmethod qfloat ((x ratio) (num-type qd-real))
 241   ;; This probably has some issues with roundoff
 242   (two-arg-/ (qfloat (numerator x) num-type)
 243              (qfloat (denominator x) num-type)))
 244
 245 (defmethod qfloat ((x ratio) (num-type qd-real))
 246   ;; This probably has some issues with roundoff
 247   (let ((top (qd-value (qfloat (numerator x) num-type)))
 248         (bot (qd-value (qfloat (denominator x) num-type))))
 249     (make-instance 'qd-real :value (div-qd top bot))))
 250
 251 #+cmu
 252 (defmethod qfloat ((x ext:double-double-float) (num-type qd-real))
 253     (make-instance 'qd-real :value (make-qd-dd x 0w0)))
 254
 255 (defmethod qfloat ((x qd-real) (num-type cl:float))
 256   (multiple-value-bind (q0 q1 q2 q3)
 257       (qd-parts (qd-value x))
 258     (cl:float (cl:+ q0 q1 q2 q3) num-type)))
 259
 260 #+cmu
 261 (defmethod qfloat ((x qd-real) (num-type ext:double-double-float))
 262   (multiple-value-bind (q0 q1 q2 q3)
 263       (qd-parts (qd-value x))
 264     (cl:+ (cl:float q0 1w0)
 265           (cl:float q1 1w0)
 266           (cl:float q2 1w0)
 267           (cl:float q3 1w0))))
 268
 269 (defmethod qfloat ((x qd-real) (num-type qd-real))
 270   x)
 271
 272 (declaim (inline float))
 273 (defun float (x num-type)
 274   (qfloat x num-type))
 275
 276 (defmethod qrealpart ((x number))
 277   (cl:realpart x))
 278 (defmethod qrealpart ((x qd-real))
 279   x)
 280 (defmethod qrealpart ((x qd-complex))
 281   (make-instance 'qd-real :value (qd-real x)))
 282 (defun realpart (x)
 283   (qrealpart x))
 284
 285 (defmethod qimagpart ((x number))
 286   (cl:imagpart x))
 287 (defmethod qimagpart ((x qd-real))
 288   (make-qd 0d0))
 289 (defmethod qimagpart ((x qd-complex))
 290   (make-instance 'qd-real :value (qd-imag x)))
 291
 292 (defun imagpart (x)
 293   (qimagpart x))
 294
 295 (defmethod qconjugate ((a number))
 296   (cl:conjugate a))
 297
 298 (defmethod qconjugate ((a qd-real))
 299   (make-instance 'qd-real :value (qd-value a)))
 300
 301 (defmethod qconjugate ((a qd-complex))
 302   (make-instance 'qd-complex
 303                  :real (qd-real a)
 304                  :imag (neg-qd (qd-imag a))))
 305
 306 (defun conjugate (z)
 307   (qconjugate z))
 308
 309 (defmethod qscale-float ((f cl:float) (n integer))
 310   (cl:scale-float f n))
 311
 312 (defmethod qscale-float ((f qd-real) (n integer))
 313   (make-instance 'qd-real :value (scale-float-qd (qd-value f) n)))
 314
 315 (declaim (inline scale-float))
 316 (defun scale-float (f n)
 317   (qscale-float f n))
 318
 319 (macrolet
 320     ((frob (op)
 321        (let ((method (intern (concatenate 'string "TWO-ARG-" (symbol-name op))))
 322              (cl-fun (find-symbol (symbol-name op) :cl))
 323              (qd-fun (intern (concatenate 'string "QD-" (symbol-name op))
 324                              (find-package :qdi))))
 325          `(progn
 326             (defmethod ,method ((a real) (b real))
 327               (,cl-fun a b))
 328             (defmethod ,method ((a qd-real) (b real))
 329               (,qd-fun (qd-value a) (make-qd-d (cl:float b 1d0))))
 330             (defmethod ,method ((a real) (b qd-real))
 331               (,qd-fun (make-qd-d (cl:float a 1d0)) (qd-value b)))
 332             (defmethod ,method ((a qd-real) (b qd-real))
 333               (,qd-fun (qd-value a) (qd-value b)))
 334             (defun ,op (number &rest more-numbers)
 335               "Returns T if its arguments are in strictly increasing order, NIL otherwise."
 336               (declare (optimize (safety 2))
 337                        (dynamic-extent more-numbers))
 338               (do* ((n number (car nlist))
 339                     (nlist more-numbers (cdr nlist)))
 340                    ((atom nlist) t)
 341                 (declare (list nlist))
 342                 (if (not (,method n (car nlist))) (return nil))))))))
 343   (frob <)
 344   (frob >)
 345   (frob <=)
 346   (frob >=))
 347
 348 (macrolet ((frob (name)
 349              (let ((method-name (intern (concatenate 'string "Q" (symbol-name name))))
 350                    (cl-name (intern (symbol-name name) :cl))
 351                    (qd-name (intern (concatenate 'string (symbol-name name) "-QD"))))
 352                `(progn
 353                  (defmethod ,method-name ((x number))
 354                    (,cl-name x))
 355                  (defmethod ,method-name ((x qd-real))
 356                    (make-instance 'qd-real :value (,qd-name (qd-value x))))
 357                  (declaim (inline ,name))
 358                  (defun ,name (x)
 359                    (,method-name x))))))
 360   (frob abs)
 361   (frob exp)
 362   (frob sin)
 363   (frob cos)
 364   (frob tan)
 365   ;;(frob asin)
 366   ;;(frob acos)
 367   (frob sinh)
 368   (frob cosh)
 369   (frob tanh)
 370   (frob asinh)
 371   ;;(frob acosh)
 372   ;;(frob atanh)
 373   )
 374
 375 (defmethod qsqrt ((x number))
 376   (cl:sqrt x))
 377
 378 (defmethod qsqrt ((x qd-real))
 379   (if (minusp x)
 380       (make-instance 'qd-complex
 381                      :real +qd-zero+
 382                      :imag (sqrt-qd (neg-qd (qd-value x))))
 383       (make-instance 'qd-real :value (sqrt-qd (qd-value x)))))
 384
 385 (defun sqrt (x)
 386   (qsqrt x))
 387
 388 (defun scalb (x n)
 389   "Compute 2^N * X without compute 2^N first (use properties of the
 390 underlying floating-point format"
 391   (declare (type qd-real x))
 392   (scale-float x n))
 393
 394 (declaim (inline qd-cssqs))
 395 (defun qd-cssqs (z)
 396   (multiple-value-bind (rho k)
 397       (qdi::hypot-aux-qd (qd-value (realpart z))
 398                          (qd-value (imagpart z)))
 399     (values (make-instance 'qd-real :value rho)
 400             k)))
 401
 402 #+nil
 403 (defmethod qabs ((z qd-complex))
 404   ;; sqrt(x^2+y^2)
 405   ;; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)
 406   (multiple-value-bind (abs^2 rho)
 407       (hypot-qd (qd-value (realpart z))
 408                 (qd-value (imagpart z)))
 409     (scale-float (make-instance 'qd-real :value (sqrt abs^2))
 410                  rho)))
 411
 412 (defmethod qabs ((z qd-complex))
 413   ;; sqrt(x^2+y^2)
 414   ;; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)
 415   (make-instance 'qd-real
 416                  :value (hypot-qd (qd-value (realpart z))
 417                                   (qd-value (imagpart z)))))
 418
 419 (defmethod qlog ((a number) &optional b)
 420   (if b
 421       (cl:log a b)
 422       (cl:log a)))
 423
 424 (defmethod qlog ((a qd-real) &optional b)
 425   (if b
 426       (/ (qlog a) (qlog b))
 427       (if (minusp (float-sign a))
 428           (make-instance 'qd-complex
 429                          :real (log-qd (abs-qd (qd-value a)))
 430                          :imag +qd-pi+)
 431           (make-instance 'qd-real :value (log-qd (qd-value a))))))
 432
 433 (declaim (inline log))
 434 (defun log (a &optional b)
 435   (qlog a b))
 436
 437
 438 (defmethod log1p ((a qd-real))
 439   (make-instance 'qd-real :value (log1p-qd (qd-value a))))
 440
 441 (defmethod qatan ((y real) &optional x)
 442   (cond (x
 443          (cond ((typep x 'qd-real)
 444                 (make-instance 'qd-real
 445                                :value (atan2-qd (qd-value y) (qd-value x))))
 446                (t
 447                 (cl:atan y x))))
 448         (t
 449          (cl:atan y))))
 450
 451 (defmethod qatan ((y qd-real) &optional x)
 452   (make-instance 'qd-real
 453                  :value
 454                  (if x
 455                      (atan2-qd (qd-value y) (qd-value x))
 456                      (atan-qd (qd-value y)))))
 457
 458 (defun atan (y &optional x)
 459   (qatan y x))
 460
 461
 462 (defmethod qexpt ((x number) (y number))
 463   (cl:expt x y))
 464
 465 (defmethod qexpt ((x qd-real) (y real))
 466   (exp (* y (log x))))
 467
 468 (defmethod qexpt ((x real) (y qd-real))
 469   (exp (* y (log x))))
 470
 471 (defmethod qexpt ((x qd-real) (y cl:complex))
 472   (exp (* (make-instance 'qd-complex
 473                          :real (qd-value (realpart y))
 474                          :imag (qd-value (imagpart y)))
 475           (log x))))
 476
 477 (defmethod qexpt ((x cl:complex) (y qd-real))
 478   (exp (* y
 479           (log (make-instance 'qd-complex
 480                               :real (qd-value (realpart x))
 481                               :imag (qd-value (imagpart x)))))))
 482
 483 (defmethod qexpt ((x qd-real) (y qd-real))
 484   ;; x^y = exp(y*log(x))
 485   (exp (* y (log x))))
 486
 487 (defmethod qexpt ((x qd-real) (y integer))
 488   (make-instance 'qd-real
 489                  :value (npow (qd-value x) y)))
 490
 491 (declaim (inline expt))
 492 (defun expt (x y)
 493   (qexpt x y))
 494
 495 \f
 496
 497 (defmethod two-arg-= ((a number) (b number))
 498   (cl:= a b))
 499 (defmethod two-arg-= ((a qd-real) (b number))
 500   (if (realp b)
 501       (qd-= (qd-value a) (make-qd-d (cl:float b 1d0)))
 502       nil))
 503 (defmethod two-arg-= ((a number) (b qd-real))
 504   (if (realp a)
 505       (qd-= (make-qd-d (cl:float a 1d0)) (qd-value b))
 506       nil))
 507
 508 (defmethod two-arg-= ((a qd-real) (b qd-real))
 509   (qd-= (qd-value a) (qd-value b)))
 510
 511 (defun = (number &rest more-numbers)
 512   "Returns T if all of its arguments are numerically equal, NIL otherwise."
 513   (declare (optimize (safety 2))
 514            (dynamic-extent more-numbers))
 515   (do ((nlist more-numbers (cdr nlist)))
 516       ((atom nlist) T)
 517     (declare (list nlist))
 518     (if (not (two-arg-= (car nlist) number))
 519         (return nil))))
 520
 521 (defun /= (number &rest more-numbers)
 522   "Returns T if no two of its arguments are numerically equal, NIL otherwise."
 523   (declare (optimize (safety 2))
 524            (dynamic-extent more-numbers))
 525   (do* ((head number (car nlist))
 526         (nlist more-numbers (cdr nlist)))
 527        ((atom nlist) t)
 528     (declare (list nlist))
 529     (unless (do* ((nl nlist (cdr nl)))
 530                  ((atom nl) T)
 531               (declare (list nl))
 532               (if (two-arg-= head (car nl))
 533                   (return nil)))
 534       (return nil))))
 535
 536 (defmethod qcomplex ((x real) &optional y)
 537   (cl:complex x (if y y 0)))
 538
 539 (defmethod qcomplex ((x qd-real) &optional y)
 540   (make-instance 'qd-complex
 541                  :real (qd-value x)
 542                  :imag (if y (qd-value y) +qd-zero+)))
 543
 544 (defun complex (x &optional (y 0))
 545   (qcomplex x y))
 546
 547 (defmethod qinteger-decode-float ((f cl:float))
 548   (cl:integer-decode-float f))
 549
 550 (defmethod qinteger-decode-float ((f qd-real))
 551   (integer-decode-qd (qd-value f)))
 552
 553 (declaim (inline integer-decode-float))
 554 (defun integer-decode-float (f)
 555   (qinteger-decode-float f))
 556
 557 (defmethod qdecode-float ((f cl:float))
 558   (cl:decode-float f))
 559
 560 (defmethod qdecode-float ((f qd-real))
 561   (multiple-value-bind (frac exp s)
 562       (decode-float-qd (qd-value f))
 563     (values (make-instance 'qd-real :value frac)
 564             exp
 565             (make-instance 'qd-real :value  s))))
 566
 567 (declaim (inline decode-float))
 568 (defun decode-float (f)
 569   (qdecode-float f))
 570
 571 (defmethod qfloor ((x real) &optional y)
 572   (if y
 573       (cl:floor x y)
 574       (cl:floor x)))
 575
 576 (defmethod qfloor ((x qd-real) &optional y)
 577   (if (and y (/= y 1))
 578       (let ((f (qfloor (/ x y))))
 579         (values f
 580                 (- x (* f y))))
 581       (let ((f (ffloor-qd (qd-value x))))
 582         (multiple-value-bind (int exp sign)
 583             (integer-decode-qd f)
 584           (values (ash (* sign int) exp)
 585                   (make-instance 'qd-real
 586                                  :value (qd-value
 587                                          (- x (make-instance 'qd-real
 588                                                              :value f)))))))))
 589
 590 (defun floor (x &optional y)
 591   (qfloor x y))
 592
 593 (defmethod qffloor ((x real) &optional y)
 594   (if y
 595       (cl:ffloor x y)
 596       (cl:ffloor x)))
 597
 598 (defmethod qffloor ((x qd-real) &optional y)
 599   (if (and y (/= y 1))
 600       (let ((f (qffloor (/ x y))))
 601         (values f
 602                 (- x (* f y))))
 603       (let ((f (make-instance 'qd-real :value (ffloor-qd (qd-value x)))))
 604         (values f
 605                 (- x f)))))
 606
 607 (defun ffloor (x &optional y)
 608   (qffloor x y))
 609
 610 (defun ceiling (x &optional y)
 611   (multiple-value-bind (f rem)
 612       (floor x y)
 613     (if (zerop rem)
 614         (values (+ f 1)
 615                 rem)
 616         (values (+ f 1)
 617                 (- rem 1)))))
 618
 619 (defun fceiling (x &optional y)
 620   (multiple-value-bind (f rem)
 621       (ffloor x y)
 622     (if (zerop rem)
 623         (values (+ f 1)
 624                 rem)
 625         (values (+ f 1)
 626                 (- rem 1)))))
 627
 628 (defun truncate (x &optional (y 1))
 629   (if (minusp x)
 630       (ceiling x y)
 631       (floor x y)))
 632
 633 (defun ftruncate (x &optional (y 1))
 634   (if (minusp x)
 635       (fceiling x y)
 636       (ffloor x y)))
 637
 638 (defmethod %unary-round ((x real))
 639   (cl::round x))
 640
 641 (defmethod %unary-round ((number qd-real))
 642   (multiple-value-bind (bits exp)
 643       (integer-decode-float number)
 644     (let* ((shifted (ash bits exp))
 645            (rounded (if (and (minusp exp)
 646                              (oddp shifted)
 647                              (not (zerop (logand bits
 648                                                  (ash 1 (- -1 exp))))))
 649                         (1+ shifted)
 650                         shifted)))
 651       (if (minusp number)
 652           (- rounded)
 653           rounded))))
 654
 655 (defun round (number &optional (divisor 1))
 656   (if (eql divisor 1)
 657       (let ((r (%unary-round number)))
 658         (values r
 659                 (- number r)))
 660       (multiple-value-bind (tru rem)
 661           (truncate number divisor)
 662         (if (zerop rem)
 663             (values tru rem)
 664             (let ((thresh (/ (abs divisor) 2)))
 665               (cond ((or (> rem thresh)
 666                          (and (= rem thresh) (oddp tru)))
 667                      (if (minusp divisor)
 668                          (values (- tru 1) (+ rem divisor))
 669                          (values (+ tru 1) (- rem divisor))))
 670                     ((let ((-thresh (- thresh)))
 671                        (or (< rem -thresh)
 672                            (and (= rem -thresh) (oddp tru))))
 673                      (if (minusp divisor)
 674                          (values (+ tru 1) (- rem divisor))
 675                          (values (- tru 1) (+ rem divisor))))
 676                     (t (values tru rem))))))))
 677
 678 (defun fround (number &optional (divisor 1))
 679   "Same as ROUND, but returns first value as a float."
 680   (multiple-value-bind (res rem)
 681       (round number divisor)
 682     (values (float res (if (floatp rem) rem 1.0)) rem)))
 683
 684 (defmethod qfloat-sign ((a real) &optional (f (float 1 a)))
 685   (cl:float-sign a f))
 686
 687 (defmethod qfloat-sign ((a qd-real) &optional f)
 688   (if f
 689       (make-instance 'qd-real
 690                      :value (mul-qd-d (abs-qd (qd-value f))
 691                                       (cl:float-sign (qd-0 (qd-value a)))))
 692       (make-instance 'qd-real :value (make-qd-d (cl:float-sign (qd-0 (qd-value a)))))))
 693
 694 (declaim (inline float-sign))
 695 (defun float-sign (n &optional float2)
 696   (qfloat-sign n float2))
 697
 698 (defun max (number &rest more-numbers)
 699   "Returns the greatest of its arguments."
 700   (declare (optimize (safety 2)) (type (or real qd-real) number)
 701            (dynamic-extent more-numbers))
 702   (dolist (real more-numbers)
 703     (when (> real number)
 704       (setq number real)))
 705   number)
 706
 707 (defun min (number &rest more-numbers)
 708   "Returns the least of its arguments."
 709   (declare (optimize (safety 2)) (type (or real qd-real) number)
 710            (dynamic-extent more-numbers))
 711   (do ((nlist more-numbers (cdr nlist))
 712        (result (the (or real qd-real) number)))
 713       ((null nlist) (return result))
 714     (declare (list nlist))
 715     (if (< (car nlist) result)
 716         (setq result (car nlist)))))
 717 \f
 718 (defmethod qasin ((x number))
 719   (cl:asin x))
 720
 721 (defmethod qasin ((x qd-real))
 722   (if (<= -1 x 1)
 723       (make-instance 'qd-real :value (asin-qd (qd-value x)))
 724       (qd-complex-asin x)))
 725
 726 (declaim (inline asin))
 727 (defun asin (x)
 728   (qasin x))
 729
 730 (defmethod qacos ((x number))
 731   (cl:acos x))
 732
 733 (defmethod qacos ((x qd-real))
 734   (cond ((> (abs x) 1)
 735          (qd-complex-acos x))
 736         (t
 737          (make-instance 'qd-real :value (acos-qd (qd-value x))))))
 738
 739 (declaim (inline acos))
 740 (defun acos (x)
 741   (qacos x))
 742
 743 (defmethod qacosh ((x number))
 744   (cl:acosh x))
 745
 746 (defmethod qacosh ((x qd-real))
 747   (if (< x 1)
 748       (qd-complex-acosh x)
 749       (make-instance 'qd-real :value (acosh-qd (qd-value x)))))
 750
 751
 752 (declaim (inline acosh))
 753 (defun acosh (x)
 754   (qacosh x))
 755
 756 (defmethod qatanh ((x number))
 757   (cl:atanh x))
 758
 759 (defmethod qatanh ((x qd-real))
 760   (if (> (abs x) 1)
 761       (qd-complex-atanh x)
 762       (make-instance 'qd-real :value (atanh-qd (qd-value x)))))
 763
 764
 765 (declaim (inline atanh))
 766 (defun atanh (x)
 767   (qatanh x))
 768
 769 (defmethod qcis ((x real))
 770   (cl:cis x))
 771
 772 (defmethod qcis ((x qd-real))
 773   (multiple-value-bind (s c)
 774       (sincos-qd (qd-value x))
 775     (make-instance 'qd-complex
 776                    :real c
 777                    :imag s)))
 778
 779 (declaim (inline cis))
 780 (defun cis (x)
 781   (qcis x))
 782
 783 (defmethod qphase ((x number))
 784   (cl:phase x))
 785
 786 (defmethod qphase ((x qd-real))
 787   (if (minusp x)
 788       (- +pi+)
 789       (make-instance 'qd-real :value (make-qd-d 0d0))))
 790
 791 (declaim (inline phase))
 792 (defun phase (x)
 793   (qphase x))
 794
 795 (defun signum (number)
 796   "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
 797   (if (zerop number)
 798       number
 799       (if (rationalp number)
 800           (if (plusp number) 1 -1)
 801           (/ number (abs number)))))
 802
 803 (defmethod random ((x cl:real) &optional (state *random-state*))
 804   (cl:random x state))
 805
 806 (defmethod random ((x qd-real) &optional (state *random-state*))
 807   (* x (make-instance 'qd-real
 808                       :value (qdi:random-qd state))))
 809 \f
 810 (define-compiler-macro + (&whole form &rest args)
 811   (if (null args)
 812       0
 813       (do ((args (cdr args) (cdr args))
 814            (res (car args)
 815                 `(two-arg-+ ,res ,(car args))))
 816           ((null args) res))))
 817
 818 (define-compiler-macro - (&whole form number &rest more-numbers)
 819   (if more-numbers
 820       (do ((nlist more-numbers (cdr nlist))
 821            (result number))
 822           ((atom nlist) result)
 823          (declare (list nlist))
 824          (setq result `(two-arg-- ,result ,(car nlist))))
 825       `(unary-minus ,number)))
 826
 827 (define-compiler-macro * (&whole form &rest args)
 828   (if (null args)
 829       1
 830       (do ((args (cdr args) (cdr args))
 831            (res (car args)
 832                 `(two-arg-* ,res ,(car args))))
 833           ((null args) res))))
 834
 835 (define-compiler-macro / (number &rest more-numbers)
 836   (if more-numbers
 837       (do ((nlist more-numbers (cdr nlist))
 838            (result number))
 839           ((atom nlist) result)
 840          (declare (list nlist))
 841          (setq result `(two-arg-/ ,result ,(car nlist))))
 842       `(unary-divide ,number)))
 843
 844 ;; Compiler macros to convert <, >, <=, and >= into multiple calls of
 845 ;; the corresponding two-arg-<foo> function.
 846 (macrolet
 847     ((frob (op)
 848        (let ((method (intern (concatenate 'string "TWO-ARG-" (symbol-name op)))))
 849          `(define-compiler-macro ,op (number &rest more-numbers)
 850             (do* ((n number (car nlist))
 851                   (nlist more-numbers (cdr nlist))
 852                   (res nil))
 853                  ((atom nlist)
 854                   `(and ,@(nreverse res)))
 855               (push `(,',method ,n ,(car nlist)) res))))))
 856   (frob <)
 857   (frob >)
 858   (frob <=)
 859   (frob >=))
 860
 861 (define-compiler-macro /= (&whole form number &rest more-numbers)
 862   ;; Convert (/= x y) to (not (two-arg-= x y)).  Should we try to
 863   ;; handle a few more cases?
 864   (if (cdr more-numbers)
 865       form
 866       `(not (two-arg-= ,number ,(car more-numbers)))))
 867
 868 \f
 869 (defun read-qd-real-or-complex (stream)
 870   (let ((c (peek-char t stream)))
 871     (cond ((char= c #\()
 872            ;; Read a QD complex
 873            (read-char stream)           ; Skip the paren
 874            (let ((real (read stream t nil t))
 875                  (imag (read stream t nil t)))
 876              (unless (char= (peek-char t stream) #\))
 877                (error "Illegal QD-COMPLEX number format"))
 878              ;; Read closing paren
 879              (read-char stream)
 880              (make-instance 'qd-complex
 881                             :real (qd-value (float real +qd-real-one+))
 882                             :imag (qd-value (float imag +qd-real-one+)))))
 883           (t
 884            (make-instance 'qd-real :value (read-qd stream))))))
 885
 886 (defun qd-class-reader (stream subchar arg)
 887   (declare (ignore subchar))
 888   (when arg
 889     (warn "Numeric argument ignored in #~DQ" arg))
 890   (read-qd-real-or-complex stream))
 891
 892 ;; Yow!  We redefine the #q reader that is in qd-io.lisp to read in
 893 ;; and make a real qd-real float, instead of the hackish
 894 ;; %qd-real.
 895 (set-dispatch-macro-character #\# #\Q #'qd-class-reader)
 896