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