Commit | Line | Data |
---|---|---|
5113f2c2 | 1 | ;;; float.el --- obsolete floating point arithmetic package |
17cd3083 | 2 | |
f2e3589a | 3 | ;; Copyright (C) 1986, 2001, 2002, 2003, 2004, 2005, |
2f043267 | 4 | ;; 2006, 2007, 2008 Free Software Foundation, Inc. |
17cd3083 JB |
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 | |
5a9dffec | 14 | ;; the Free Software Foundation; either version 3, or (at your option) |
17cd3083 JB |
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 | |
3a35cf56 LK |
24 | ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, |
25 | ;; Boston, MA 02110-1301, USA. | |
17cd3083 JB |
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: | |
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 | (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, | |
235 | nil 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 | |
249 | second, 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, | |
254 | nil 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 | |
259 | second, 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, | |
264 | nil 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, | |
316 | like 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. | |
327 | Optional 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. | |
396 | Accepts a decimal string in scientific notation, with exponent preceded | |
397 | by either E or e. Only the six most significant digits of the integer | |
398 | and fractional parts are used; only the first two digits of the exponent | |
399 | are used. Negative signs preceding both the decimal number and the exponent | |
400 | are 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 | |
2b8f53de | 427 | (f (* (+ (string-to-number digit-string) |
17cd3083 JB |
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 | ||
2b8f53de | 437 | (setq expt (+ (* (string-to-number |
17cd3083 JB |
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 | ||
cbee283d | 459 | ;; arch-tag: cc0c89c6-5718-49af-978e-585f6b14e347 |
17cd3083 | 460 | ;;; float.el ends here |