1 ;;; float.el --- floating point arithmetic package.
3 ;; Author: Bill Rosenblatt
5 ;; Last-Modified: 16 Mar 1992
7 ;; Copyright (C) 1986 Free Software Foundation, Inc.
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is free software; you can redistribute it and/or modify
12 ;; it under the terms of the GNU General Public License as published by
13 ;; the Free Software Foundation; either version 2, or (at your option)
16 ;; GNU Emacs is distributed in the hope that it will be useful,
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ;; GNU General Public License for more details.
21 ;; You should have received a copy of the GNU General Public License
22 ;; along with GNU Emacs; see the file COPYING. If not, write to
23 ;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
27 ;; Floating point numbers are represented by dot-pairs (mant . exp)
28 ;; where mant is the 24-bit signed integral mantissa and exp is the
31 ;; Emacs LISP supports a 24-bit signed integer data type, which has a
32 ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
33 ;; This gives six significant decimal digit accuracy. Exponents can
34 ;; be anything in the range -(2**23) to +(2**23)-1.
37 ;; function f converts from integer to floating point
38 ;; function string-to-float converts from string to floating point
39 ;; function fint converts a floating point to integer (with truncation)
40 ;; function float-to-string converts from floating point to string
43 ;; - Exponents outside of the range of +/-100 or so will cause certain
44 ;; functions (especially conversion routines) to take forever.
45 ;; - Very little checking is done for fixed point overflow/underflow.
46 ;; - No checking is done for over/underflow of the exponent
47 ;; (hardly necessary when exponent can be 2**23).
56 ;; fundamental implementation constants
58 "Base of exponent in this floating point representation.")
60 (defconst mantissa-bits
24
61 "Number of significant bits in this floating point representation.")
63 (defconst decimal-digits
6
64 "Number of decimal digits expected to be accurate.")
66 (defconst expt-digits
2
67 "Maximum permitted digits in a scientific notation exponent.")
70 (defconst maxbit
(1- mantissa-bits
)
71 "Number of highest bit")
73 (defconst mantissa-maxval
(1- (ash 1 maxbit
))
74 "Maximum permissable value of mantissa")
76 (defconst mantissa-minval
(ash 1 maxbit
)
77 "Minimum permissable value of mantissa")
79 (defconst floating-point-regexp
80 "^[ \t]*\\(-?\\)\\([0-9]*\\)\
81 \\(\\.\\([0-9]*\\)\\|\\)\
82 \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
83 "Regular expression to match floating point numbers. Extract matches:
87 8 - minus sign for power of ten
91 (defconst high-bit-mask
(ash 1 maxbit
)
92 "Masks all bits except the high-order (sign) bit.")
94 (defconst second-bit-mask
(ash 1 (1- maxbit
))
95 "Masks all bits except the highest-order magnitude bit")
97 ;; various useful floating point constants
100 (setq _f1
/2 '(4194304 . -
23))
102 (setq _f1
'(4194304 . -
22))
104 (setq _f10
'(5242880 . -
19))
106 ;; support for decimal conversion routines
107 (setq powers-of-10
(make-vector (1+ decimal-digits
) _f1
))
108 (aset powers-of-10
1 _f10
)
109 (aset powers-of-10
2 '(6553600 . -
16))
110 (aset powers-of-10
3 '(8192000 . -
13))
111 (aset powers-of-10
4 '(5120000 . -
9))
112 (aset powers-of-10
5 '(6400000 . -
6))
113 (aset powers-of-10
6 '(8000000 . -
3))
115 (setq all-decimal-digs-minval
(aref powers-of-10
(1- decimal-digits
))
116 highest-power-of-10
(aref powers-of-10 decimal-digits
))
118 (defun fashl (fnum) ; floating-point arithmetic shift left
119 (cons (ash (car fnum
) 1) (1- (cdr fnum
))))
121 (defun fashr (fnum) ; floating point arithmetic shift right
122 (cons (ash (car fnum
) -
1) (1+ (cdr fnum
))))
124 (defun normalize (fnum)
125 (if (> (car fnum
) 0) ; make sure next-to-highest bit is set
126 (while (zerop (logand (car fnum
) second-bit-mask
))
127 (setq fnum
(fashl fnum
)))
128 (if (< (car fnum
) 0) ; make sure highest bit is set
129 (while (zerop (logand (car fnum
) high-bit-mask
))
130 (setq fnum
(fashl fnum
)))
131 (setq fnum _f0
))) ; "standard 0"
134 (defun abs (n) ; integer absolute value
135 (if (>= n
0) n
(- n
)))
137 (defun fabs (fnum) ; re-normalize after taking abs value
138 (normalize (cons (abs (car fnum
)) (cdr fnum
))))
140 (defun xor (a b
) ; logical exclusive or
141 (and (or a b
) (not (and a b
))))
143 (defun same-sign (a b
) ; two f-p numbers have same sign?
144 (not (xor (natnump (car a
)) (natnump (car b
)))))
146 (defun extract-match (str i
) ; used after string-match
148 (substring str
(match-beginning i
) (match-end i
))
151 ;; support for the multiplication function
152 (setq halfword-bits
(/ mantissa-bits
2) ; bits in a halfword
153 masklo
(1- (ash 1 halfword-bits
)) ; isolate the lower halfword
154 maskhi
(lognot masklo
) ; isolate the upper halfword
155 round-limit
(ash 1 (/ halfword-bits
2)))
157 (defun hihalf (n) ; return high halfword, shifted down
158 (ash (logand n maskhi
) (- halfword-bits
)))
160 (defun lohalf (n) ; return low halfword
165 ;; Arithmetic functions
167 "Returns the sum of two floating point numbers."
168 (let ((f1 (fmax a1 a2
))
170 (if (same-sign a1 a2
)
171 (setq f1
(fashr f1
) ; shift right to avoid overflow
174 (cons (+ (car f1
) (ash (car f2
) (- (cdr f2
) (cdr f1
))))
177 (defun f- (a1 &optional a2
) ; unary or binary minus
178 "Returns the difference of two floating point numbers."
181 (normalize (cons (- (car a1
)) (cdr a1
)))))
183 (defun f* (a1 a2
) ; multiply in halfword chunks
184 "Returns the product of two floating point numbers."
185 (let* ((i1 (car (fabs a1
)))
187 (sign (not (same-sign a1 a2
)))
188 (prodlo (+ (hihalf (* (lohalf i1
) (lohalf i2
)))
189 (lohalf (* (hihalf i1
) (lohalf i2
)))
190 (lohalf (* (lohalf i1
) (hihalf i2
)))))
191 (prodhi (+ (* (hihalf i1
) (hihalf i2
))
192 (hihalf (* (hihalf i1
) (lohalf i2
)))
193 (hihalf (* (lohalf i1
) (hihalf i2
)))
195 (if (> (lohalf prodlo
) round-limit
)
196 (setq prodhi
(1+ prodhi
))) ; round off truncated bits
198 (cons (if sign
(- prodhi
) prodhi
)
199 (+ (cdr (fabs a1
)) (cdr (fabs a2
)) mantissa-bits
)))))
201 (defun f/ (a1 a2
) ; SLOW subtract-and-shift algorithm
202 "Returns the quotient of two floating point numbers."
203 (if (zerop (car a2
)) ; if divide by 0
204 (signal 'arith-error
(list "attempt to divide by zero" a1 a2
))
205 (let ((bits (1- maxbit
))
207 (dividend (car (fabs a1
)))
208 (divisor (car (fabs a2
)))
209 (sign (not (same-sign a1 a2
))))
210 (while (natnump bits
)
211 (if (< (- dividend divisor
) 0)
212 (setq quotient
(ash quotient
1))
213 (setq quotient
(1+ (ash quotient
1))
214 dividend
(- dividend divisor
)))
215 (setq dividend
(ash dividend
1)
218 (cons (if sign
(- quotient
) quotient
)
219 (- (cdr (fabs a1
)) (cdr (fabs a2
)) (1- maxbit
)))))))
222 "Returns the remainder of first floating point number divided by second."
223 (f- a1
(f* (ftrunc (f/ a1 a2
)) a2
)))
226 ;; Comparison functions
228 "Returns t if two floating point numbers are equal, nil otherwise."
232 "Returns t if first floating point number is greater than second,
234 (cond ((and (natnump (car a1
)) (< (car a2
) 0))
235 t
) ; a1 nonnegative, a2 negative
236 ((and (> (car a1
) 0) (<= (car a2
) 0))
237 t
) ; a1 positive, a2 nonpositive
238 ((and (<= (car a1
) 0) (natnump (car a2
)))
239 nil
) ; a1 nonpos, a2 nonneg
240 ((/= (cdr a1
) (cdr a2
)) ; same signs. exponents differ
241 (> (cdr a1
) (cdr a2
))) ; compare the mantissas.
243 (> (car a1
) (car a2
))))) ; same exponents.
246 "Returns t if first floating point number is greater than or equal to
247 second, nil otherwise."
248 (or (f> a1 a2
) (f= a1 a2
)))
251 "Returns t if first floating point number is less than second,
256 "Returns t if first floating point number is less than or equal to
257 second, nil otherwise."
261 "Returns t if first floating point number is not equal to second,
266 "Returns the minimum of two floating point numbers."
267 (if (f< a1 a2
) a1 a2
))
270 "Returns the maximum of two floating point numbers."
271 (if (f> a1 a2
) a1 a2
))
274 "Returns t if the floating point number is zero, nil otherwise."
278 "Returns t if the arg is a floating point number, nil otherwise."
279 (and (consp fnum
) (integerp (car fnum
)) (integerp (cdr fnum
))))
281 ;; Conversion routines
283 "Convert the integer argument to floating point, like a C cast operator."
284 (normalize (cons int
'0)))
286 (defun int-to-hex-string (int)
287 "Convert the integer argument to a C-style hexadecimal string."
290 (hex-chars "0123456789ABCDEF"))
291 (while (<= shiftval
0)
292 (setq str
(concat str
(char-to-string
294 (logand (lsh int shiftval
) 15))))
295 shiftval
(+ shiftval
4)))
298 (defun ftrunc (fnum) ; truncate fractional part
299 "Truncate the fractional part of a floating point number."
300 (cond ((natnump (cdr fnum
)) ; it's all integer, return number as is
302 ((<= (cdr fnum
) (- maxbit
)) ; it's all fractional, return 0
304 (t ; otherwise mask out fractional bits
305 (let ((mant (car fnum
)) (exp (cdr fnum
)))
307 (cons (if (natnump mant
) ; if negative, use absolute value
308 (ash (ash mant exp
) (- exp
))
309 (- (ash (ash (- mant
) exp
) (- exp
))))
312 (defun fint (fnum) ; truncate and convert to integer
313 "Convert the floating point number to integer, with truncation,
314 like a C cast operator."
315 (let* ((tf (ftrunc fnum
)) (tint (car tf
)) (texp (cdr tf
)))
316 (cond ((>= texp mantissa-bits
) ; too high, return "maxint"
318 ((<= texp
(- mantissa-bits
)) ; too low, return "minint"
321 (ash tint texp
))))) ; shift so that exponent is 0
323 (defun float-to-string (fnum &optional sci
)
324 "Convert the floating point number to a decimal string.
325 Optional second argument non-nil means use scientific notation."
326 (let* ((value (fabs fnum
)) (sign (< (car fnum
) 0))
327 (power 0) (result 0) (str "")
328 (temp 0) (pow10 _f1
))
332 (if (f>= value _f1
) ; find largest power of 10 <= value
333 (progn ; value >= 1, power is positive
334 (while (f<= (setq temp
(f* pow10 highest-power-of-10
)) value
)
336 power
(+ power decimal-digits
)))
337 (while (f<= (setq temp
(f* pow10 _f10
)) value
)
340 (progn ; value < 1, power is negative
341 (while (f> (setq temp
(f/ pow10 highest-power-of-10
)) value
)
343 power
(- power decimal-digits
)))
344 (while (f> pow10 value
)
345 (setq pow10
(f/ pow10 _f10
)
347 ; get value in range 100000 to 999999
348 (setq value
(f* (f/ value pow10
) all-decimal-digs-minval
)
349 result
(ftrunc value
))
351 (if (f> (f- value result
) _f1
/2) ; round up if remainder > 0.5
352 (setq int
(1+ (fint result
)))
353 (setq int
(fint result
)))
354 (setq str
(int-to-string int
))
356 (setq power
(1+ power
))))
358 (if sci
; scientific notation
359 (setq str
(concat (substring str
0 1) "." (substring str
1)
360 "E" (int-to-string power
)))
362 ; regular decimal string
363 (cond ((>= power
(1- decimal-digits
))
364 ; large power, append zeroes
365 (let ((zeroes (- power decimal-digits
)))
366 (while (natnump zeroes
)
367 (setq str
(concat str
"0")
368 zeroes
(1- zeroes
)))))
370 ; negative power, prepend decimal
371 ((< power
0) ; point and zeroes
372 (let ((zeroes (- (- power
) 2)))
373 (while (natnump zeroes
)
374 (setq str
(concat "0" str
)
376 (setq str
(concat "0." str
))))
378 (t ; in range, insert decimal point
380 (substring str
0 (1+ power
))
382 (substring str
(1+ power
)))))))
384 (if sign
; if negative, prepend minus sign
389 ;; string to float conversion.
390 ;; accepts scientific notation, but ignores anything after the first two
391 ;; digits of the exponent.
392 (defun string-to-float (str)
393 "Convert the string to a floating point number.
394 Accepts a decimal string in scientific notation, with exponent preceded
395 by either E or e. Only the six most significant digits of the integer
396 and fractional parts are used; only the first two digits of the exponent
397 are used. Negative signs preceding both the decimal number and the exponent
400 (if (string-match floating-point-regexp str
0)
403 ; calculate the mantissa
404 (let* ((int-subst (extract-match str
2))
405 (fract-subst (extract-match str
4))
406 (digit-string (concat int-subst fract-subst
))
407 (mant-sign (equal (extract-match str
1) "-"))
408 (leading-0s 0) (round-up nil
))
410 ; get rid of leading 0's
411 (setq power
(- (length int-subst
) decimal-digits
))
412 (while (and (< leading-0s
(length digit-string
))
413 (= (aref digit-string leading-0s
) ?
0))
414 (setq leading-0s
(1+ leading-0s
)))
415 (setq power
(- power leading-0s
)
416 digit-string
(substring digit-string leading-0s
))
418 ; if more than 6 digits, round off
419 (if (> (length digit-string
) decimal-digits
)
420 (setq round-up
(>= (aref digit-string decimal-digits
) ?
5)
421 digit-string
(substring digit-string
0 decimal-digits
))
422 (setq power
(+ power
(- decimal-digits
(length digit-string
)))))
424 ; round up and add minus sign, if necessary
425 (f (* (+ (string-to-int digit-string
)
427 (if mant-sign -
1 1))))
429 ; calculate the exponent (power of ten)
430 (let* ((expt-subst (extract-match str
9))
431 (expt-sign (equal (extract-match str
8) "-"))
432 (expt 0) (chunks 0) (tens 0) (exponent _f1
)
435 (setq expt
(+ (* (string-to-int
436 (substring expt-subst
0
437 (min expt-digits
(length expt-subst
))))
440 (if (< expt
0) ; if power of 10 negative
441 (setq expt
(- expt
) ; take abs val of exponent
442 func
'f
/)) ; and set up to divide, not multiply
444 (setq chunks
(/ expt decimal-digits
)
445 tens
(% expt decimal-digits
))
446 ; divide or multiply by "chunks" of 10**6
448 (setq exponent
(funcall func exponent highest-power-of-10
)
450 ; divide or multiply by remaining power of ten
451 (funcall func exponent
(aref powers-of-10 tens
)))))
453 _f0
)) ; if invalid, return 0
457 ;;; float.el ends here