1;;; float.el --- obsolete floating point arithmetic package 2 3;; Copyright (C) 1986, 2001, 2002, 2003, 2004, 2005, 4;; 2006, 2007 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 2, 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: 861 - minus sign 872 - integer part 884 - fractional part 898 - minus sign for power of ten 909 - 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, 235nil 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 249second, 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, 254nil 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 259second, 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, 264nil 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, 316like 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. 327Optional 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. 396Accepts a decimal string in scientific notation, with exponent preceded 397by either E or e. Only the six most significant digits of the integer 398and fractional parts are used; only the first two digits of the exponent 399are used. Negative signs preceding both the decimal number and the exponent 400are 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 461