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