lisp/obsolete/float.el

   1 ;;; float.el --- obsolete floating point arithmetic package
   2
   3 ;; Copyright (C) 1986, 2001, 2002, 2003, 2004, 2005,
   4 ;;   2006, 2007, 2008 Free Software Foundation, Inc.
   5
   6 ;; Author: Bill Rosenblatt
   7 ;; Maintainer: FSF
   8 ;; Keywords: extensions
   9
  10 ;; This file is part of GNU Emacs.
  11
  12 ;; GNU Emacs is free software; you can redistribute it and/or modify
  13 ;; it under the terms of the GNU General Public License as published by
  14 ;; the Free Software Foundation; either version 3, or (at your option)
  15 ;; any later version.
  16
  17 ;; GNU Emacs is distributed in the hope that it will be useful,
  18 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20 ;; GNU General Public License for more details.
  21
  22 ;; You should have received a copy of the GNU General Public License
  23 ;; along with GNU Emacs; see the file COPYING.  If not, write to the
  24 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  25 ;; Boston, MA 02110-1301, USA.
  26
  27 ;;; Commentary:
  28
  29 ;; Floating point numbers are represented by dot-pairs (mant . exp)
  30 ;; where mant is the 24-bit signed integral mantissa and exp is the
  31 ;; base 2 exponent.
  32 ;;
  33 ;; Emacs LISP supports a 24-bit signed integer data type, which has a
  34 ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
  35 ;; This gives six significant decimal digit accuracy.  Exponents can
  36 ;; be anything in the range -(2**23) to +(2**23)-1.
  37 ;;
  38 ;; User interface:
  39 ;; function f converts from integer to floating point
  40 ;; function string-to-float converts from string to floating point
  41 ;; function fint converts a floating point to integer (with truncation)
  42 ;; function float-to-string converts from floating point to string
  43 ;;
  44 ;; Caveats:
  45 ;; -  Exponents outside of the range of +/-100 or so will cause certain
  46 ;;    functions (especially conversion routines) to take forever.
  47 ;; -  Very little checking is done for fixed point overflow/underflow.
  48 ;; -  No checking is done for over/underflow of the exponent
  49 ;;    (hardly necessary when exponent can be 2**23).
  50 ;;
  51 ;;
  52 ;; Bill Rosenblatt
  53 ;; June 20, 1986
  54 ;;
  55
  56 ;;; Code:
  57
  58 ;; fundamental implementation constants
  59 (defconst exp-base 2
  60   "Base of exponent in this floating point representation.")
  61
  62 (defconst mantissa-bits 24
  63   "Number of significant bits in this floating point representation.")
  64
  65 (defconst decimal-digits 6
  66   "Number of decimal digits expected to be accurate.")
  67
  68 (defconst expt-digits 2
  69   "Maximum permitted digits in a scientific notation exponent.")
  70
  71 ;; other constants
  72 (defconst maxbit (1- mantissa-bits)
  73   "Number of highest bit")
  74
  75 (defconst mantissa-maxval (1- (ash 1 maxbit))
  76   "Maximum permissible value of mantissa")
  77
  78 (defconst mantissa-minval (ash 1 maxbit)
  79   "Minimum permissible value of mantissa")
  80
  81 (defconst floating-point-regexp
  82   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
  83 \\(\\.\\([0-9]*\\)\\|\\)\
  84 \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  85   "Regular expression to match floating point numbers.  Extract matches:
  86 1 - minus sign
  87 2 - integer part
  88 4 - fractional part
  89 8 - minus sign for power of ten
  90 9 - power of ten
  91 ")
  92
  93 (defconst high-bit-mask (ash 1 maxbit)
  94   "Masks all bits except the high-order (sign) bit.")
  95
  96 (defconst second-bit-mask (ash 1 (1- maxbit))
  97   "Masks all bits except the highest-order magnitude bit")
  98
  99 ;; various useful floating point constants
 100 (defconst _f0 '(0 . 1))
 101
 102 (defconst _f1/2 '(4194304 . -23))
 103
 104 (defconst _f1 '(4194304 . -22))
 105
 106 (defconst _f10 '(5242880 . -19))
 107
 108 ;; support for decimal conversion routines
 109 (defvar powers-of-10 (make-vector (1+ decimal-digits) _f1))
 110 (aset powers-of-10 1 _f10)
 111 (aset powers-of-10 2 '(6553600 . -16))
 112 (aset powers-of-10 3 '(8192000 . -13))
 113 (aset powers-of-10 4 '(5120000 . -9))
 114 (aset powers-of-10 5 '(6400000 . -6))
 115 (aset powers-of-10 6 '(8000000 . -3))
 116
 117 (defconst all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits)))
 118 (defconst highest-power-of-10 (aref powers-of-10 decimal-digits))
 119
 120 (defun fashl (fnum)                     ; floating-point arithmetic shift left
 121   (cons (ash (car fnum) 1) (1- (cdr fnum))))
 122
 123 (defun fashr (fnum)                     ; floating point arithmetic shift right
 124   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
 125
 126 (defun normalize (fnum)
 127   (if (> (car fnum) 0)                  ; make sure next-to-highest bit is set
 128       (while (zerop (logand (car fnum) second-bit-mask))
 129         (setq fnum (fashl fnum)))
 130     (if (< (car fnum) 0)                ; make sure highest bit is set
 131         (while (zerop (logand (car fnum) high-bit-mask))
 132           (setq fnum (fashl fnum)))
 133       (setq fnum _f0)))                 ; "standard 0"
 134   fnum)
 135
 136 (defun abs (n)                          ; integer absolute value
 137   (if (>= n 0) n (- n)))
 138
 139 (defun fabs (fnum)                      ; re-normalize after taking abs value
 140   (normalize (cons (abs (car fnum)) (cdr fnum))))
 141
 142 (defun xor (a b)                        ; logical exclusive or
 143   (and (or a b) (not (and a b))))
 144
 145 (defun same-sign (a b)                  ; two f-p numbers have same sign?
 146   (not (xor (natnump (car a)) (natnump (car b)))))
 147
 148 (defun extract-match (str i)            ; used after string-match
 149   (condition-case ()
 150       (substring str (match-beginning i) (match-end i))
 151     (error "")))
 152
 153 ;; support for the multiplication function
 154 (defconst halfword-bits (/ mantissa-bits 2)) ; bits in a halfword
 155 (defconst masklo (1- (ash 1 halfword-bits))) ; isolate the lower halfword
 156 (defconst maskhi (lognot masklo))       ; isolate the upper halfword
 157 (defconst round-limit (ash 1 (/ halfword-bits 2)))
 158
 159 (defun hihalf (n)                       ; return high halfword, shifted down
 160   (ash (logand n maskhi) (- halfword-bits)))
 161
 162 (defun lohalf (n)                       ; return low halfword
 163   (logand n masklo))
 164
 165 ;; Visible functions
 166
 167 ;; Arithmetic functions
 168 (defun f+ (a1 a2)
 169   "Returns the sum of two floating point numbers."
 170   (let ((f1 (fmax a1 a2))
 171         (f2 (fmin a1 a2)))
 172     (if (same-sign a1 a2)
 173         (setq f1 (fashr f1)             ; shift right to avoid overflow
 174               f2 (fashr f2)))
 175     (normalize
 176      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
 177            (cdr f1)))))
 178
 179 (defun f- (a1 &optional a2)             ; unary or binary minus
 180   "Returns the difference of two floating point numbers."
 181   (if a2
 182       (f+ a1 (f- a2))
 183     (normalize (cons (- (car a1)) (cdr a1)))))
 184
 185 (defun f* (a1 a2)                       ; multiply in halfword chunks
 186   "Returns the product of two floating point numbers."
 187   (let* ((i1 (car (fabs a1)))
 188          (i2 (car (fabs a2)))
 189          (sign (not (same-sign a1 a2)))
 190          (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
 191                     (lohalf (* (hihalf i1) (lohalf i2)))
 192                     (lohalf (* (lohalf i1) (hihalf i2)))))
 193          (prodhi (+ (* (hihalf i1) (hihalf i2))
 194                     (hihalf (* (hihalf i1) (lohalf i2)))
 195                     (hihalf (* (lohalf i1) (hihalf i2)))
 196                     (hihalf prodlo))))
 197     (if (> (lohalf prodlo) round-limit)
 198         (setq prodhi (1+ prodhi)))      ; round off truncated bits
 199     (normalize
 200      (cons (if sign (- prodhi) prodhi)
 201            (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
 202
 203 (defun f/ (a1 a2)                       ; SLOW subtract-and-shift algorithm
 204   "Returns the quotient of two floating point numbers."
 205   (if (zerop (car a2))                  ; if divide by 0
 206       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
 207     (let ((bits (1- maxbit))
 208           (quotient 0)
 209           (dividend (car (fabs a1)))
 210           (divisor (car (fabs a2)))
 211           (sign (not (same-sign a1 a2))))
 212       (while (natnump bits)
 213         (if (< (- dividend divisor) 0)
 214             (setq quotient (ash quotient 1))
 215           (setq quotient (1+ (ash quotient 1))
 216                 dividend (- dividend divisor)))
 217         (setq dividend (ash dividend 1)
 218               bits (1- bits)))
 219       (normalize
 220        (cons (if sign (- quotient) quotient)
 221              (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
 222
 223 (defun f% (a1 a2)
 224   "Returns the remainder of first floating point number divided by second."
 225   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
 226
 227
 228 ;; Comparison functions
 229 (defun f= (a1 a2)
 230   "Returns t if two floating point numbers are equal, nil otherwise."
 231   (equal a1 a2))
 232
 233 (defun f> (a1 a2)
 234   "Returns t if first floating point number is greater than second,
 235 nil otherwise."
 236   (cond ((and (natnump (car a1)) (< (car a2) 0))
 237          t)                             ; a1 nonnegative, a2 negative
 238         ((and (> (car a1) 0) (<= (car a2) 0))
 239          t)                             ; a1 positive, a2 nonpositive
 240         ((and (<= (car a1) 0) (natnump (car a2)))
 241          nil)                           ; a1 nonpos, a2 nonneg
 242         ((/= (cdr a1) (cdr a2))         ; same signs.  exponents differ
 243          (> (cdr a1) (cdr a2)))         ; compare the mantissas.
 244         (t
 245          (> (car a1) (car a2)))))       ; same exponents.
 246
 247 (defun f>= (a1 a2)
 248   "Returns t if first floating point number is greater than or equal to
 249 second, nil otherwise."
 250   (or (f> a1 a2) (f= a1 a2)))
 251
 252 (defun f< (a1 a2)
 253   "Returns t if first floating point number is less than second,
 254 nil otherwise."
 255   (not (f>= a1 a2)))
 256
 257 (defun f<= (a1 a2)
 258   "Returns t if first floating point number is less than or equal to
 259 second, nil otherwise."
 260   (not (f> a1 a2)))
 261
 262 (defun f/= (a1 a2)
 263   "Returns t if first floating point number is not equal to second,
 264 nil otherwise."
 265   (not (f= a1 a2)))
 266
 267 (defun fmin (a1 a2)
 268   "Returns the minimum of two floating point numbers."
 269   (if (f< a1 a2) a1 a2))
 270
 271 (defun fmax (a1 a2)
 272   "Returns the maximum of two floating point numbers."
 273   (if (f> a1 a2) a1 a2))
 274
 275 (defun fzerop (fnum)
 276   "Returns t if the floating point number is zero, nil otherwise."
 277   (= (car fnum) 0))
 278
 279 (defun floatp (fnum)
 280   "Returns t if the arg is a floating point number, nil otherwise."
 281   (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
 282
 283 ;; Conversion routines
 284 (defun f (int)
 285   "Convert the integer argument to floating point, like a C cast operator."
 286   (normalize (cons int '0)))
 287
 288 (defun int-to-hex-string (int)
 289   "Convert the integer argument to a C-style hexadecimal string."
 290   (let ((shiftval -20)
 291         (str "0x")
 292         (hex-chars "0123456789ABCDEF"))
 293     (while (<= shiftval 0)
 294       (setq str (concat str (char-to-string
 295                         (aref hex-chars
 296                               (logand (lsh int shiftval) 15))))
 297             shiftval (+ shiftval 4)))
 298     str))
 299
 300 (defun ftrunc (fnum)                    ; truncate fractional part
 301   "Truncate the fractional part of a floating point number."
 302   (cond ((natnump (cdr fnum))           ; it's all integer, return number as is
 303          fnum)
 304         ((<= (cdr fnum) (- maxbit))     ; it's all fractional, return 0
 305          '(0 . 1))
 306         (t                              ; otherwise mask out fractional bits
 307          (let ((mant (car fnum)) (exp (cdr fnum)))
 308            (normalize
 309             (cons (if (natnump mant)    ; if negative, use absolute value
 310                       (ash (ash mant exp) (- exp))
 311                     (- (ash (ash (- mant) exp) (- exp))))
 312                   exp))))))
 313
 314 (defun fint (fnum)                      ; truncate and convert to integer
 315   "Convert the floating point number to integer, with truncation,
 316 like a C cast operator."
 317   (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
 318     (cond ((>= texp mantissa-bits)      ; too high, return "maxint"
 319            mantissa-maxval)
 320           ((<= texp (- mantissa-bits))  ; too low, return "minint"
 321            mantissa-minval)
 322           (t                            ; in range
 323            (ash tint texp)))))          ; shift so that exponent is 0
 324
 325 (defun float-to-string (fnum &optional sci)
 326   "Convert the floating point number to a decimal string.
 327 Optional second argument non-nil means use scientific notation."
 328   (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
 329          (power 0) (result 0) (str "")
 330          (temp 0) (pow10 _f1))
 331
 332     (if (f= fnum _f0)
 333         "0"
 334       (if (f>= value _f1)                       ; find largest power of 10 <= value
 335           (progn                                ; value >= 1, power is positive
 336             (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
 337               (setq pow10 temp
 338                     power (+ power decimal-digits)))
 339             (while (f<= (setq temp (f* pow10 _f10)) value)
 340               (setq pow10 temp
 341                     power (1+ power))))
 342         (progn                          ; value < 1, power is negative
 343           (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
 344             (setq pow10 temp
 345                   power (- power decimal-digits)))
 346           (while (f> pow10 value)
 347             (setq pow10 (f/ pow10 _f10)
 348                   power (1- power)))))
 349                                           ; get value in range 100000 to 999999
 350       (setq value (f* (f/ value pow10) all-decimal-digs-minval)
 351             result (ftrunc value))
 352       (let (int)
 353         (if (f> (f- value result) _f1/2)        ; round up if remainder > 0.5
 354             (setq int (1+ (fint result)))
 355           (setq int (fint result)))
 356         (setq str (int-to-string int))
 357         (if (>= int 1000000)
 358             (setq power (1+ power))))
 359
 360       (if sci                           ; scientific notation
 361           (setq str (concat (substring str 0 1) "." (substring str 1)
 362                             "E" (int-to-string power)))
 363
 364                                           ; regular decimal string
 365         (cond ((>= power (1- decimal-digits))
 366                                           ; large power, append zeroes
 367                (let ((zeroes (- power decimal-digits)))
 368                  (while (natnump zeroes)
 369                    (setq str (concat str "0")
 370                          zeroes (1- zeroes)))))
 371
 372                                           ; negative power, prepend decimal
 373               ((< power 0)              ; point and zeroes
 374                (let ((zeroes (- (- power) 2)))
 375                  (while (natnump zeroes)
 376                    (setq str (concat "0" str)
 377                          zeroes (1- zeroes)))
 378                  (setq str (concat "0." str))))
 379
 380               (t                                ; in range, insert decimal point
 381                (setq str (concat
 382                           (substring str 0 (1+ power))
 383                           "."
 384                           (substring str (1+ power)))))))
 385
 386       (if sign                          ; if negative, prepend minus sign
 387           (concat "-" str)
 388         str))))
 389
 390
 391 ;; string to float conversion.
 392 ;; accepts scientific notation, but ignores anything after the first two
 393 ;; digits of the exponent.
 394 (defun string-to-float (str)
 395   "Convert the string to a floating point number.
 396 Accepts a decimal string in scientific notation, with exponent preceded
 397 by either E or e.  Only the six most significant digits of the integer
 398 and fractional parts are used; only the first two digits of the exponent
 399 are used.  Negative signs preceding both the decimal number and the exponent
 400 are recognized."
 401
 402   (if (string-match floating-point-regexp str 0)
 403       (let (power)
 404         (f*
 405          ; calculate the mantissa
 406          (let* ((int-subst (extract-match str 2))
 407                 (fract-subst (extract-match str 4))
 408                 (digit-string (concat int-subst fract-subst))
 409                 (mant-sign (equal (extract-match str 1) "-"))
 410                 (leading-0s 0) (round-up nil))
 411
 412            ; get rid of leading 0's
 413            (setq power (- (length int-subst) decimal-digits))
 414            (while (and (< leading-0s (length digit-string))
 415                        (= (aref digit-string leading-0s) ?0))
 416              (setq leading-0s (1+ leading-0s)))
 417            (setq power (- power leading-0s)
 418                  digit-string (substring digit-string leading-0s))
 419
 420            ; if more than 6 digits, round off
 421            (if (> (length digit-string) decimal-digits)
 422                (setq round-up (>= (aref digit-string decimal-digits) ?5)
 423                      digit-string (substring digit-string 0 decimal-digits))
 424              (setq power (+ power (- decimal-digits (length digit-string)))))
 425
 426            ; round up and add minus sign, if necessary
 427            (f (* (+ (string-to-number digit-string)
 428                     (if round-up 1 0))
 429                  (if mant-sign -1 1))))
 430
 431          ; calculate the exponent (power of ten)
 432          (let* ((expt-subst (extract-match str 9))
 433                 (expt-sign (equal (extract-match str 8) "-"))
 434                 (expt 0) (chunks 0) (tens 0) (exponent _f1)
 435                 (func 'f*))
 436
 437            (setq expt (+ (* (string-to-number
 438                              (substring expt-subst 0
 439                                         (min expt-digits (length expt-subst))))
 440                             (if expt-sign -1 1))
 441                          power))
 442            (if (< expt 0)               ; if power of 10 negative
 443                (setq expt (- expt)      ; take abs val of exponent
 444                      func 'f/))         ; and set up to divide, not multiply
 445
 446            (setq chunks (/ expt decimal-digits)
 447                  tens (% expt decimal-digits))
 448            ; divide or multiply by "chunks" of 10**6
 449            (while (> chunks 0)
 450              (setq exponent (funcall func exponent highest-power-of-10)
 451                    chunks (1- chunks)))
 452            ; divide or multiply by remaining power of ten
 453            (funcall func exponent (aref powers-of-10 tens)))))
 454
 455     _f0))                               ; if invalid, return 0
 456
 457 (provide 'float)
 458
 459 ;; arch-tag: cc0c89c6-5718-49af-978e-585f6b14e347
 460 ;;; float.el ends here