1 ;; Copyright (C) 1986 Free Software Foundation, Inc.
2 ;; Author Bill Rosenblatt
4 ;; This file is part of GNU Emacs.
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)
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.
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.
20 ;; Floating point arithmetic package.
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
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.
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
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).
49 ;; fundamental implementation constants
51 "Base of exponent in this floating point representation.")
53 (defconst mantissa-bits 24
54 "Number of significant bits in this floating point representation.")
56 (defconst decimal-digits 6
57 "Number of decimal digits expected to be accurate.")
59 (defconst expt-digits 2
60 "Maximum permitted digits in a scientific notation exponent.")
63 (defconst maxbit (1- mantissa-bits)
64 "Number of highest bit")
66 (defconst mantissa-maxval (1- (ash 1 maxbit))
67 "Maximum permissable value of mantissa")
69 (defconst mantissa-minval (ash 1 maxbit)
70 "Minimum permissable value of mantissa")
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:
80 8 - minus sign for power of ten
84 (defconst high-bit-mask (ash 1 maxbit)
85 "Masks all bits except the high-order (sign) bit.")
87 (defconst second-bit-mask (ash 1 (1- maxbit))
88 "Masks all bits except the highest-order magnitude bit")
90 ;; various useful floating point constants
93 (setq _f1/2 '(4194304 . -23))
95 (setq _f1 '(4194304 . -22))
97 (setq _f10 '(5242880 . -19))
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))
108 (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
109 highest-power-of-10 (aref powers-of-10 decimal-digits))
111 (defun fashl (fnum) ; floating-point arithmetic shift left
112 (cons (ash (car fnum) 1) (1- (cdr fnum))))
114 (defun fashr (fnum) ; floating point arithmetic shift right
115 (cons (ash (car fnum) -1) (1+ (cdr fnum))))
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"
127 (defun abs (n) ; integer absolute value
128 (if (>= n 0) n (- n)))
130 (defun fabs (fnum) ; re-normalize after taking abs value
131 (normalize (cons (abs (car fnum)) (cdr fnum))))
133 (defun xor (a b) ; logical exclusive or
134 (and (or a b) (not (and a b))))
136 (defun same-sign (a b) ; two f-p numbers have same sign?
137 (not (xor (natnump (car a)) (natnump (car b)))))
139 (defun extract-match (str i) ; used after string-match
141 (substring str (match-beginning i) (match-end i))
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)))
150 (defun hihalf (n) ; return high halfword, shifted down
151 (ash (logand n maskhi) (- halfword-bits)))
153 (defun lohalf (n) ; return low halfword
158 ;; Arithmetic functions
160 "Returns the sum of two floating point numbers."
161 (let ((f1 (fmax a1 a2))
163 (if (same-sign a1 a2)
164 (setq f1 (fashr f1) ; shift right to avoid overflow
167 (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
170 (defun f- (a1 &optional a2) ; unary or binary minus
171 "Returns the difference of two floating point numbers."
174 (normalize (cons (- (car a1)) (cdr a1)))))
176 (defun f* (a1 a2) ; multiply in halfword chunks
177 "Returns the product of two floating point numbers."
178 (let* ((i1 (car (fabs a1)))
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)))
188 (if (> (lohalf prodlo) round-limit)
189 (setq prodhi (1+ prodhi))) ; round off truncated bits
191 (cons (if sign (- prodhi) prodhi)
192 (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
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))
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)
211 (cons (if sign (- quotient) quotient)
212 (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
215 "Returns the remainder of first floating point number divided by second."
216 (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
219 ;; Comparison functions
221 "Returns t if two floating point numbers are equal, nil otherwise."
225 "Returns t if first floating point number is greater than second,
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.
236 (> (car a1) (car a2))))) ; same exponents.
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)))
244 "Returns t if first floating point number is less than second,
249 "Returns t if first floating point number is less than or equal to
250 second, nil otherwise."
254 "Returns t if first floating point number is not equal to second,
259 "Returns the minimum of two floating point numbers."
260 (if (f< a1 a2) a1 a2))
263 "Returns the maximum of two floating point numbers."
264 (if (f> a1 a2) a1 a2))
267 "Returns t if the floating point number is zero, nil otherwise."
271 "Returns t if the arg is a floating point number, nil otherwise."
272 (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
274 ;; Conversion routines
276 "Convert the integer argument to floating point, like a C cast operator."
277 (normalize (cons int '0)))
279 (defun int-to-hex-string (int)
280 "Convert the integer argument to a C-style hexadecimal string."
283 (hex-chars "0123456789ABCDEF"))
284 (while (<= shiftval 0)
285 (setq str (concat str (char-to-string
287 (logand (lsh int shiftval) 15))))
288 shiftval (+ shiftval 4)))
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
295 ((<= (cdr fnum) (- maxbit)) ; it's all fractional, return 0
297 (t ; otherwise mask out fractional bits
298 (let ((mant (car fnum)) (exp (cdr fnum)))
300 (cons (if (natnump mant) ; if negative, use absolute value
301 (ash (ash mant exp) (- exp))
302 (- (ash (ash (- mant) exp) (- exp))))
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"
311 ((<= texp (- mantissa-bits)) ; too low, return "minint"
314 (ash tint texp))))) ; shift so that exponent is 0
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))
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)
329 power (+ power decimal-digits)))
330 (while (f<= (setq temp (f* pow10 _f10)) value)
333 (progn ; value < 1, power is negative
334 (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
336 power (- power decimal-digits)))
337 (while (f> pow10 value)
338 (setq pow10 (f/ pow10 _f10)
340 ; get value in range 100000 to 999999
341 (setq value (f* (f/ value pow10) all-decimal-digs-minval)
342 result (ftrunc value))
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))
349 (setq power (1+ power))))
351 (if sci ; scientific notation
352 (setq str (concat (substring str 0 1) "." (substring str 1)
353 "E" (int-to-string power)))
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)))))
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)
369 (setq str (concat "0." str))))
371 (t ; in range, insert decimal point
373 (substring str 0 (1+ power))
375 (substring str (1+ power)))))))
377 (if sign ; if negative, prepend minus sign
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.
387 Accepts a decimal string in scientific notation,
388 with exponent preceded by either E or e.
389 Only the 6 most significant digits of the integer and fractional parts
390 are used; only the first two digits of the exponent are used.
391 Negative signs preceding both the decimal number and the exponent
394 (if (string-match floating-point-regexp str 0)
397 ; calculate the mantissa
398 (let* ((int-subst (extract-match str 2))
399 (fract-subst (extract-match str 4))
400 (digit-string (concat int-subst fract-subst))
401 (mant-sign (equal (extract-match str 1) "-"))
402 (leading-0s 0) (round-up nil))
404 ; get rid of leading 0's
405 (setq power (- (length int-subst) decimal-digits))
406 (while (and (< leading-0s (length digit-string))
407 (= (aref digit-string leading-0s) ?0))
408 (setq leading-0s (1+ leading-0s)))
409 (setq power (- power leading-0s)
410 digit-string (substring digit-string leading-0s))
412 ; if more than 6 digits, round off
413 (if (> (length digit-string) decimal-digits)
414 (setq round-up (>= (aref digit-string decimal-digits) ?5)
415 digit-string (substring digit-string 0 decimal-digits))
416 (setq power (+ power (- decimal-digits (length digit-string)))))
418 ; round up and add minus sign, if necessary
419 (f (* (+ (string-to-int digit-string)
421 (if mant-sign -1 1))))
423 ; calculate the exponent (power of ten)
424 (let* ((expt-subst (extract-match str 9))
425 (expt-sign (equal (extract-match str 8) "-"))
426 (expt 0) (chunks 0) (tens 0) (exponent _f1)
429 (setq expt (+ (* (string-to-int
430 (substring expt-subst 0
431 (min expt-digits (length expt-subst))))
434 (if (< expt 0) ; if power of 10 negative
435 (setq expt (- expt) ; take abs val of exponent
436 func 'f/)) ; and set up to divide, not multiply
438 (setq chunks (/ expt decimal-digits)
439 tens (% expt decimal-digits))
440 ; divide or multiply by "chunks" of 10**6
442 (setq exponent (funcall func exponent highest-power-of-10)
444 ; divide or multiply by remaining power of ten
445 (funcall func exponent (aref powers-of-10 tens)))))
447 _f0)) ; if invalid, return 0