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