numlib.lsp

   1 ;; (c) Copyright Taiichi Yuasa and Masami Hagiya, 1984.  All rights reserved.
   2 ;; Copying of this file is authorized to users who have executed the true and
   3 ;; proper "License Agreement for Kyoto Common LISP" with SIGLISP.
   4
   5 ;;;;    numlib.lsp
   6 ;;;;
   7 ;;;;                           number routines
   8
   9
  10 (in-package 'lisp)
  11 (export
  12  '(isqrt abs phase signum cis asin acos sinh cosh tanh
  13    asinh acosh atanh
  14    rational rationalize
  15    ffloor fround ftruncate fceiling
  16    lognand lognor logandc1 logandc2 logorc1 logorc2
  17    lognot logtest
  18    byte byte-size byte-position
  19    ldb ldb-test mask-field dpb deposit-field
  20    ))
  21
  22
  23 (in-package 'system)
  24
  25
  26 (proclaim '(optimize (safety 2) (space 3)))
  27
  28
  29 (defconstant imag-one #C(0.0f0 1.0f0))
  30
  31
  32 (defun isqrt (i)
  33        (unless (and (integerp i) (>= i 0))
  34                (error "~S is not a non-negative integer." i))
  35        (if (zerop i)
  36            0
  37            (let ((n (integer-length i)))
  38                 (do ((x (ash 1 (ceiling n 2)))
  39                      (y))
  40                     (nil)
  41                     (setq y (floor i x))
  42                     (when (<= x y)
  43                           (return x))
  44                     (setq x (floor (+ x y) 2))))))
  45
  46 (defun abs (x)
  47        (if (complexp x)
  48            (sqrt (+ (* (realpart x) (realpart x))
  49                     (* (imagpart x) (imagpart x))))
  50            (if (minusp x)
  51                (- x)
  52                x)))
  53
  54 (defun phase (x)
  55        (atan (imagpart x) (realpart x)))
  56
  57 (defun signum (x) (if (zerop x) x (/ x (abs x))))
  58
  59 (defun cis (x) (exp (* imag-one x)))
  60
  61
  62 #|
  63 (defun asin (x)
  64        (let ((c (- (* imag-one
  65                       (log (+ (* imag-one x)
  66                               (sqrt (- 1.0f0 (* x x)))))))))
  67             (if (or (not (complexp x)) (zerop (imagpart c)))
  68                 (realpart c)
  69                 c)))
  70
  71 (defun acos (x)
  72        (let ((c (- (* imag-one
  73                       (log (+ x (* imag-one
  74                                    (sqrt (- 1.0f0 (* x x))))))))))
  75             (if (or (not (complexp x)) (zerop (imagpart c)))
  76                 (realpart c)
  77                 c)))
  78 |#
  79
  80 (defun sinh (x) (/ (- (exp x) (exp (- x))) 2.0f0))
  81 (defun cosh (x) (/ (+ (exp x) (exp (- x))) 2.0f0))
  82 (defun tanh (x) (/ (sinh x) (cosh x)))
  83
  84 (defun asinh (x) (log (+ x (sqrt (+ 1.0f0 (* x x))))))
  85 (defun acosh (x)
  86        (log (+ x
  87                (* (1+ x)
  88                     (sqrt (/ (1- x) (1+ x)))))))
  89 (defun atanh (x)
  90        (when (or (= x 1.0f0) (= x -1.0f0))
  91              (error "The argument, ~s, is a logarithmic singularity.~
  92                     ~%Don't be foolish, GLS."
  93                     x))
  94        (log (/ (1+ x) (sqrt (- 1.0f0 (* x x))))))
  95
  96
  97 (defun rational (x)
  98   (etypecase x
  99     (float
 100       (multiple-value-bind (i e s) (integer-decode-float x)
 101                            (if (>= s 0)
 102                                (* i (expt (float-radix x) e))
 103                              (- (* i (expt (float-radix x) e))))))
 104     (rational x)))
 105
 106
 107 (setf (symbol-function 'rationalize) (symbol-function 'rational))
 108
 109 ;; although the following is correct code in that it approximates the
 110 ;; x to within eps, it does not preserve (eql (float (rationalize x) x) x)
 111 ;; since the test for eql is more strict than the float-epsilon
 112
 113 ;;; Rationalize originally by Skef Wholey.
 114 ;;; Obtained from Daniel L. Weinreb.
 115 ;(defun rationalize (x)
 116 ;  (typecase x
 117 ;    (rational x)
 118 ;    (short-float (rationalize-float x short-float-epsilon 1.0s0))
 119 ;    (long-float (rationalize-float x long-float-epsilon 1.0d0))
 120 ;    (otherwise (error "~S is neither rational nor float." x))))
 121 ;
 122 ;(defun rationalize-float (x eps one)
 123 ;  (cond ((minusp x) (- (rationalize (- x))))
 124 ;        ((zerop x) 0)
 125 ;        (t (let ((y ())
 126 ;                 (a ()))
 127 ;             (do ((xx x (setq y (/ one
 128 ;                                   (- xx (float a x)))))
 129 ;                  (num (setq a (truncate x))
 130 ;                       (+ (* (setq a (truncate y)) num) onum))
 131 ;                  (den 1 (+ (* a den) oden))
 132 ;                  (onum 1 num)
 133 ;                  (oden 0 den))
 134 ;                 ((and (not (zerop den))
 135 ;                       (not (> (abs (/ (- x (/ (float num x)
 136 ;                                               (float den x)))
 137 ;                                       x))
 138 ;                               eps)))
 139 ;                  (/ num den)))))))
 140
 141
 142 (defun ffloor (x &optional (y 1.0f0))
 143        (multiple-value-bind (i r) (floor (float x) (float y))
 144         (values (float i r) r)))
 145
 146 (defun fceiling (x &optional (y 1.0f0))
 147        (multiple-value-bind (i r) (ceiling (float x) (float y))
 148         (values (float i r) r)))
 149
 150 (defun ftruncate (x &optional (y 1.0f0))
 151        (multiple-value-bind (i r) (truncate (float x) (float y))
 152         (values (float i r) r)))
 153
 154 (defun fround (x &optional (y 1.0f0))
 155        (multiple-value-bind (i r) (round (float x) (float y))
 156         (values (float i r) r)))
 157
 158
 159 (defun lognand (x y) (boole boole-nand x y))
 160 (defun lognor (x y) (boole boole-nor x y))
 161 (defun logandc1 (x y) (boole boole-andc1 x y))
 162 (defun logandc2 (x y) (boole boole-andc2 x y))
 163 (defun logorc1 (x y) (boole boole-orc1 x y))
 164 (defun logorc2 (x y) (boole boole-orc2 x y))
 165
 166 (defun lognot (x) (logxor -1 x))
 167 (defun logtest (x y) (not (zerop (logand x y))))
 168
 169
 170 (defun byte (size position)
 171   (cons size position))
 172
 173 (defun byte-size (bytespec)
 174   (car bytespec))
 175
 176 (defun byte-position (bytespec)
 177   (cdr bytespec))
 178
 179 (defun ldb (bytespec integer)
 180   (logandc2 (ash integer (- (byte-position bytespec)))
 181             (- (ash 1 (byte-size bytespec)))))
 182
 183 (defun ldb-test (bytespec integer)
 184   (not (zerop (ldb bytespec integer))))
 185
 186 (defun mask-field (bytespec integer)
 187   (ash (ldb bytespec integer) (byte-position bytespec)))
 188
 189 (defun dpb (newbyte bytespec integer)
 190   (logxor integer
 191           (mask-field bytespec integer)
 192           (ash (logandc2 newbyte
 193                          (- (ash 1 (byte-size bytespec))))
 194                (byte-position bytespec))))
 195
 196 (defun deposit-field (newbyte bytespec integer)
 197   (dpb (ash newbyte (- (byte-position bytespec))) bytespec integer))