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