59
|
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
|
67
|
49 (provide 'float)
|
|
50
|
59
|
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 (defconst mantissa-minval (ash 1 maxbit)
|
|
72 "Minimum permissable value of mantissa")
|
|
73
|
|
74 (defconst floating-point-regexp
|
|
75 "^[ \t]*\\(-?\\)\\([0-9]*\\)\
|
|
76 \\(\\.\\([0-9]*\\)\\|\\)\
|
|
77 \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
|
|
78 "Regular expression to match floating point numbers. Extract matches:
|
|
79 1 - minus sign
|
|
80 2 - integer part
|
|
81 4 - fractional part
|
|
82 8 - minus sign for power of ten
|
|
83 9 - power of ten
|
|
84 ")
|
|
85
|
|
86 (defconst high-bit-mask (ash 1 maxbit)
|
|
87 "Masks all bits except the high-order (sign) bit.")
|
|
88
|
|
89 (defconst second-bit-mask (ash 1 (1- maxbit))
|
|
90 "Masks all bits except the highest-order magnitude bit")
|
|
91
|
|
92 ;; various useful floating point constants
|
|
93 (setq _f0 '(0 . 1))
|
|
94
|
|
95 (setq _f1/2 '(4194304 . -23))
|
|
96
|
|
97 (setq _f1 '(4194304 . -22))
|
|
98
|
|
99 (setq _f10 '(5242880 . -19))
|
|
100
|
|
101 ;; support for decimal conversion routines
|
|
102 (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
|
|
103 (aset powers-of-10 1 _f10)
|
|
104 (aset powers-of-10 2 '(6553600 . -16))
|
|
105 (aset powers-of-10 3 '(8192000 . -13))
|
|
106 (aset powers-of-10 4 '(5120000 . -9))
|
|
107 (aset powers-of-10 5 '(6400000 . -6))
|
|
108 (aset powers-of-10 6 '(8000000 . -3))
|
|
109
|
|
110 (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
|
|
111 highest-power-of-10 (aref powers-of-10 decimal-digits))
|
|
112
|
|
113 (defun fashl (fnum) ; floating-point arithmetic shift left
|
|
114 (cons (ash (car fnum) 1) (1- (cdr fnum))))
|
|
115
|
|
116 (defun fashr (fnum) ; floating point arithmetic shift right
|
|
117 (cons (ash (car fnum) -1) (1+ (cdr fnum))))
|
|
118
|
|
119 (defun normalize (fnum)
|
|
120 (if (> (car fnum) 0) ; make sure next-to-highest bit is set
|
|
121 (while (zerop (logand (car fnum) second-bit-mask))
|
|
122 (setq fnum (fashl fnum)))
|
|
123 (if (< (car fnum) 0) ; make sure highest bit is set
|
|
124 (while (zerop (logand (car fnum) high-bit-mask))
|
|
125 (setq fnum (fashl fnum)))
|
|
126 (setq fnum _f0))) ; "standard 0"
|
|
127 fnum)
|
|
128
|
|
129 (defun abs (n) ; integer absolute value
|
|
130 (if (>= n 0) n (- n)))
|
|
131
|
|
132 (defun fabs (fnum) ; re-normalize after taking abs value
|
|
133 (normalize (cons (abs (car fnum)) (cdr fnum))))
|
|
134
|
|
135 (defun xor (a b) ; logical exclusive or
|
|
136 (and (or a b) (not (and a b))))
|
|
137
|
|
138 (defun same-sign (a b) ; two f-p numbers have same sign?
|
|
139 (not (xor (natnump (car a)) (natnump (car b)))))
|
|
140
|
|
141 (defun extract-match (str i) ; used after string-match
|
|
142 (condition-case ()
|
|
143 (substring str (match-beginning i) (match-end i))
|
|
144 (error "")))
|
|
145
|
|
146 ;; support for the multiplication function
|
|
147 (setq halfword-bits (/ mantissa-bits 2) ; bits in a halfword
|
|
148 masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
|
|
149 maskhi (lognot masklo) ; isolate the upper halfword
|
|
150 round-limit (ash 1 (/ halfword-bits 2)))
|
|
151
|
|
152 (defun hihalf (n) ; return high halfword, shifted down
|
|
153 (ash (logand n maskhi) (- halfword-bits)))
|
|
154
|
|
155 (defun lohalf (n) ; return low halfword
|
|
156 (logand n masklo))
|
|
157
|
|
158 ;; Visible functions
|
|
159
|
|
160 ;; Arithmetic functions
|
|
161 (defun f+ (a1 a2)
|
|
162 "Returns the sum of two floating point numbers."
|
|
163 (let ((f1 (fmax a1 a2))
|
|
164 (f2 (fmin a1 a2)))
|
|
165 (if (same-sign a1 a2)
|
|
166 (setq f1 (fashr f1) ; shift right to avoid overflow
|
|
167 f2 (fashr f2)))
|
|
168 (normalize
|
|
169 (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
|
|
170 (cdr f1)))))
|
|
171
|
|
172 (defun f- (a1 &optional a2) ; unary or binary minus
|
|
173 "Returns the difference of two floating point numbers."
|
|
174 (if a2
|
|
175 (f+ a1 (f- a2))
|
|
176 (normalize (cons (- (car a1)) (cdr a1)))))
|
|
177
|
|
178 (defun f* (a1 a2) ; multiply in halfword chunks
|
|
179 "Returns the product of two floating point numbers."
|
|
180 (let* ((i1 (car (fabs a1)))
|
|
181 (i2 (car (fabs a2)))
|
|
182 (sign (not (same-sign a1 a2)))
|
|
183 (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
|
|
184 (lohalf (* (hihalf i1) (lohalf i2)))
|
|
185 (lohalf (* (lohalf i1) (hihalf i2)))))
|
|
186 (prodhi (+ (* (hihalf i1) (hihalf i2))
|
|
187 (hihalf (* (hihalf i1) (lohalf i2)))
|
|
188 (hihalf (* (lohalf i1) (hihalf i2)))
|
|
189 (hihalf prodlo))))
|
|
190 (if (> (lohalf prodlo) round-limit)
|
|
191 (setq prodhi (1+ prodhi))) ; round off truncated bits
|
|
192 (normalize
|
|
193 (cons (if sign (- prodhi) prodhi)
|
|
194 (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
|
|
195
|
|
196 (defun f/ (a1 a2) ; SLOW subtract-and-shift algorithm
|
|
197 "Returns the quotient of two floating point numbers."
|
|
198 (if (zerop (car a2)) ; if divide by 0
|
|
199 (signal 'arith-error (list "attempt to divide by zero" a1 a2))
|
|
200 (let ((bits (1- maxbit))
|
|
201 (quotient 0)
|
|
202 (dividend (car (fabs a1)))
|
|
203 (divisor (car (fabs a2)))
|
|
204 (sign (not (same-sign a1 a2))))
|
|
205 (while (natnump bits)
|
|
206 (if (< (- dividend divisor) 0)
|
|
207 (setq quotient (ash quotient 1))
|
|
208 (setq quotient (1+ (ash quotient 1))
|
|
209 dividend (- dividend divisor)))
|
|
210 (setq dividend (ash dividend 1)
|
|
211 bits (1- bits)))
|
|
212 (normalize
|
|
213 (cons (if sign (- quotient) quotient)
|
|
214 (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
|
|
215
|
|
216 (defun f% (a1 a2)
|
|
217 "Returns the remainder of first floating point number divided by second."
|
|
218 (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
|
|
219
|
|
220
|
|
221 ;; Comparison functions
|
|
222 (defun f= (a1 a2)
|
|
223 "Returns t if two floating point numbers are equal, nil otherwise."
|
|
224 (equal a1 a2))
|
|
225
|
|
226 (defun f> (a1 a2)
|
|
227 "Returns t if first floating point number is greater than second,
|
|
228 nil otherwise."
|
|
229 (cond ((and (natnump (car a1)) (< (car a2) 0))
|
|
230 t) ; a1 nonnegative, a2 negative
|
|
231 ((and (> (car a1) 0) (<= (car a2) 0))
|
|
232 t) ; a1 positive, a2 nonpositive
|
|
233 ((and (<= (car a1) 0) (natnump (car a2)))
|
|
234 nil) ; a1 nonpos, a2 nonneg
|
|
235 ((/= (cdr a1) (cdr a2)) ; same signs. exponents differ
|
|
236 (> (cdr a1) (cdr a2))) ; compare the mantissas.
|
|
237 (t
|
|
238 (> (car a1) (car a2))))) ; same exponents.
|
|
239
|
|
240 (defun f>= (a1 a2)
|
|
241 "Returns t if first floating point number is greater than or equal to
|
|
242 second, nil otherwise."
|
|
243 (or (f> a1 a2) (f= a1 a2)))
|
|
244
|
|
245 (defun f< (a1 a2)
|
|
246 "Returns t if first floating point number is less than second,
|
|
247 nil otherwise."
|
|
248 (not (f>= a1 a2)))
|
|
249
|
|
250 (defun f<= (a1 a2)
|
|
251 "Returns t if first floating point number is less than or equal to
|
|
252 second, nil otherwise."
|
|
253 (not (f> a1 a2)))
|
|
254
|
|
255 (defun f/= (a1 a2)
|
|
256 "Returns t if first floating point number is not equal to second,
|
|
257 nil otherwise."
|
|
258 (not (f= a1 a2)))
|
|
259
|
|
260 (defun fmin (a1 a2)
|
|
261 "Returns the minimum of two floating point numbers."
|
|
262 (if (f< a1 a2) a1 a2))
|
|
263
|
|
264 (defun fmax (a1 a2)
|
|
265 "Returns the maximum of two floating point numbers."
|
|
266 (if (f> a1 a2) a1 a2))
|
|
267
|
|
268 (defun fzerop (fnum)
|
|
269 "Returns t if the floating point number is zero, nil otherwise."
|
|
270 (= (car fnum) 0))
|
|
271
|
|
272 (defun floatp (fnum)
|
|
273 "Returns t if the arg is a floating point number, nil otherwise."
|
|
274 (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
|
|
275
|
|
276 ;; Conversion routines
|
|
277 (defun f (int)
|
|
278 "Convert the integer argument to floating point, like a C cast operator."
|
|
279 (normalize (cons int '0)))
|
|
280
|
|
281 (defun int-to-hex-string (int)
|
|
282 "Convert the integer argument to a C-style hexadecimal string."
|
|
283 (let ((shiftval -20)
|
|
284 (str "0x")
|
|
285 (hex-chars "0123456789ABCDEF"))
|
|
286 (while (<= shiftval 0)
|
|
287 (setq str (concat str (char-to-string
|
|
288 (aref hex-chars
|
|
289 (logand (lsh int shiftval) 15))))
|
|
290 shiftval (+ shiftval 4)))
|
|
291 str))
|
|
292
|
|
293 (defun ftrunc (fnum) ; truncate fractional part
|
|
294 "Truncate the fractional part of a floating point number."
|
|
295 (cond ((natnump (cdr fnum)) ; it's all integer, return number as is
|
|
296 fnum)
|
|
297 ((<= (cdr fnum) (- maxbit)) ; it's all fractional, return 0
|
|
298 '(0 . 1))
|
|
299 (t ; otherwise mask out fractional bits
|
|
300 (let ((mant (car fnum)) (exp (cdr fnum)))
|
|
301 (normalize
|
|
302 (cons (if (natnump mant) ; if negative, use absolute value
|
|
303 (ash (ash mant exp) (- exp))
|
|
304 (- (ash (ash (- mant) exp) (- exp))))
|
|
305 exp))))))
|
|
306
|
|
307 (defun fint (fnum) ; truncate and convert to integer
|
|
308 "Convert the floating point number to integer, with truncation,
|
|
309 like a C cast operator."
|
|
310 (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
|
|
311 (cond ((>= texp mantissa-bits) ; too high, return "maxint"
|
|
312 mantissa-maxval)
|
|
313 ((<= texp (- mantissa-bits)) ; too low, return "minint"
|
|
314 mantissa-minval)
|
|
315 (t ; in range
|
|
316 (ash tint texp))))) ; shift so that exponent is 0
|
|
317
|
|
318 (defun float-to-string (fnum &optional sci)
|
|
319 "Convert the floating point number to a decimal string.
|
|
320 Optional second argument non-nil means use scientific notation."
|
|
321 (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
|
|
322 (power 0) (result 0) (str "")
|
|
323 (temp 0) (pow10 _f1))
|
|
324
|
|
325 (if (f= fnum _f0)
|
|
326 "0"
|
|
327 (if (f>= value _f1) ; find largest power of 10 <= value
|
|
328 (progn ; value >= 1, power is positive
|
|
329 (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
|
|
330 (setq pow10 temp
|
|
331 power (+ power decimal-digits)))
|
|
332 (while (f<= (setq temp (f* pow10 _f10)) value)
|
|
333 (setq pow10 temp
|
|
334 power (1+ power))))
|
|
335 (progn ; value < 1, power is negative
|
|
336 (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
|
|
337 (setq pow10 temp
|
|
338 power (- power decimal-digits)))
|
|
339 (while (f> pow10 value)
|
|
340 (setq pow10 (f/ pow10 _f10)
|
|
341 power (1- power)))))
|
|
342 ; get value in range 100000 to 999999
|
|
343 (setq value (f* (f/ value pow10) all-decimal-digs-minval)
|
|
344 result (ftrunc value))
|
|
345 (let (int)
|
|
346 (if (f> (f- value result) _f1/2) ; round up if remainder > 0.5
|
|
347 (setq int (1+ (fint result)))
|
|
348 (setq int (fint result)))
|
|
349 (setq str (int-to-string int))
|
|
350 (if (>= int 1000000)
|
|
351 (setq power (1+ power))))
|
|
352
|
|
353 (if sci ; scientific notation
|
|
354 (setq str (concat (substring str 0 1) "." (substring str 1)
|
|
355 "E" (int-to-string power)))
|
|
356
|
|
357 ; regular decimal string
|
|
358 (cond ((>= power (1- decimal-digits))
|
|
359 ; large power, append zeroes
|
|
360 (let ((zeroes (- power decimal-digits)))
|
|
361 (while (natnump zeroes)
|
|
362 (setq str (concat str "0")
|
|
363 zeroes (1- zeroes)))))
|
|
364
|
|
365 ; negative power, prepend decimal
|
|
366 ((< power 0) ; point and zeroes
|
|
367 (let ((zeroes (- (- power) 2)))
|
|
368 (while (natnump zeroes)
|
|
369 (setq str (concat "0" str)
|
|
370 zeroes (1- zeroes)))
|
|
371 (setq str (concat "0." str))))
|
|
372
|
|
373 (t ; in range, insert decimal point
|
|
374 (setq str (concat
|
|
375 (substring str 0 (1+ power))
|
|
376 "."
|
|
377 (substring str (1+ power)))))))
|
|
378
|
|
379 (if sign ; if negative, prepend minus sign
|
|
380 (concat "-" str)
|
|
381 str))))
|
|
382
|
|
383
|
|
384 ;; string to float conversion.
|
|
385 ;; accepts scientific notation, but ignores anything after the first two
|
|
386 ;; digits of the exponent.
|
|
387 (defun string-to-float (str)
|
|
388 "Convert the string to a floating point number.
|
|
389 Accepts a decimal string in scientific notation,
|
|
390 with exponent preceded by either E or e.
|
|
391 Only the 6 most significant digits of the integer and fractional parts
|
|
392 are used; only the first two digits of the exponent are used.
|
|
393 Negative signs preceding both the decimal number and the exponent
|
|
394 are recognized."
|
|
395
|
|
396 (if (string-match floating-point-regexp str 0)
|
|
397 (let (power)
|
|
398 (f*
|
|
399 ; calculate the mantissa
|
|
400 (let* ((int-subst (extract-match str 2))
|
|
401 (fract-subst (extract-match str 4))
|
|
402 (digit-string (concat int-subst fract-subst))
|
|
403 (mant-sign (equal (extract-match str 1) "-"))
|
|
404 (leading-0s 0) (round-up nil))
|
|
405
|
|
406 ; get rid of leading 0's
|
|
407 (setq power (- (length int-subst) decimal-digits))
|
|
408 (while (and (< leading-0s (length digit-string))
|
|
409 (= (aref digit-string leading-0s) ?0))
|
|
410 (setq leading-0s (1+ leading-0s)))
|
|
411 (setq power (- power leading-0s)
|
|
412 digit-string (substring digit-string leading-0s))
|
|
413
|
|
414 ; if more than 6 digits, round off
|
|
415 (if (> (length digit-string) decimal-digits)
|
|
416 (setq round-up (>= (aref digit-string decimal-digits) ?5)
|
|
417 digit-string (substring digit-string 0 decimal-digits))
|
|
418 (setq power (+ power (- decimal-digits (length digit-string)))))
|
|
419
|
|
420 ; round up and add minus sign, if necessary
|
|
421 (f (* (+ (string-to-int digit-string)
|
|
422 (if round-up 1 0))
|
|
423 (if mant-sign -1 1))))
|
|
424
|
|
425 ; calculate the exponent (power of ten)
|
|
426 (let* ((expt-subst (extract-match str 9))
|
|
427 (expt-sign (equal (extract-match str 8) "-"))
|
|
428 (expt 0) (chunks 0) (tens 0) (exponent _f1)
|
|
429 (func 'f*))
|
|
430
|
|
431 (setq expt (+ (* (string-to-int
|
|
432 (substring expt-subst 0
|
|
433 (min expt-digits (length expt-subst))))
|
|
434 (if expt-sign -1 1))
|
|
435 power))
|
|
436 (if (< expt 0) ; if power of 10 negative
|
|
437 (setq expt (- expt) ; take abs val of exponent
|
|
438 func 'f/)) ; and set up to divide, not multiply
|
|
439
|
|
440 (setq chunks (/ expt decimal-digits)
|
|
441 tens (% expt decimal-digits))
|
|
442 ; divide or multiply by "chunks" of 10**6
|
|
443 (while (> chunks 0)
|
|
444 (setq exponent (funcall func exponent highest-power-of-10)
|
|
445 chunks (1- chunks)))
|
|
446 ; divide or multiply by remaining power of ten
|
|
447 (funcall func exponent (aref powers-of-10 tens)))))
|
|
448
|
|
449 _f0)) ; if invalid, return 0
|
|
450
|
|
451
|