![[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 [CSRG BSD Unix] / 43BSDReno / contrib / emacs-18.55 / lisp|
Return to float.el CVS log | Up to [CSRG BSD Unix] / 43BSDReno / contrib / emacs-18.55 / 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 distributed in the hope that it will be useful,
7: ;; but WITHOUT ANY WARRANTY. No author or distributor
8: ;; accepts responsibility to anyone for the consequences of using it
9: ;; or for whether it serves any particular purpose or works at all,
10: ;; unless he says so in writing. Refer to the GNU Emacs General Public
11: ;; License for full details.
12:
13: ;; Everyone is granted permission to copy, modify and redistribute
14: ;; GNU Emacs, but only under the conditions described in the
15: ;; GNU Emacs General Public License. A copy of this license is
16: ;; supposed to have been given to you along with GNU Emacs so you
17: ;; can know your rights and responsibilities. It should be in a
18: ;; file named COPYING. Among other things, the copyright notice
19: ;; and this notice must be preserved on all copies.
20:
21: ;; Floating point arithmetic package.
22: ;;
23: ;; Floating point numbers are represented by dot-pairs (mant . exp)
24: ;; where mant is the 24-bit signed integral mantissa and exp is the
25: ;; base 2 exponent.
26: ;;
27: ;; Emacs LISP supports a 24-bit signed integer data type, which has a
28: ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
29: ;; This gives six significant decimal digit accuracy. Exponents can
30: ;; be anything in the range -(2**23) to +(2**23)-1.
31: ;;
32: ;; User interface:
33: ;; function f converts from integer to floating point
34: ;; function string-to-float converts from string to floating point
35: ;; function fint converts a floating point to integer (with truncation)
36: ;; function float-to-string converts from floating point to string
37: ;;
38: ;; Caveats:
39: ;; - Exponents outside of the range of +/-100 or so will cause certain
40: ;; functions (especially conversion routines) to take forever.
41: ;; - Very little checking is done for fixed point overflow/underflow.
42: ;; - No checking is done for over/underflow of the exponent
43: ;; (hardly necessary when exponent can be 2**23).
44: ;;
45: ;;
46: ;; Bill Rosenblatt
47: ;; June 20, 1986
48: ;;
49:
50: ;; fundamental implementation constants
51: (defconst exp-base 2
52: "Base of exponent in this floating point representation.")
53:
54: (defconst mantissa-bits 24
55: "Number of significant bits in this floating point representation.")
56:
57: (defconst decimal-digits 6
58: "Number of decimal digits expected to be accurate.")
59:
60: (defconst expt-digits 2
61: "Maximum permitted digits in a scientific notation exponent.")
62:
63: ;; other constants
64: (defconst maxbit (1- mantissa-bits)
65: "Number of highest bit")
66:
67: (defconst mantissa-maxval (1- (ash 1 maxbit))
68: "Maximum permissable value of mantissa")
69:
70: (defconst mantissa-minval (ash 1 maxbit)
71: "Minimum permissable value of mantissa")
72:
73: (defconst floating-point-regexp
74: "^[ \t]*\\(-?\\)\\([0-9]*\\)\
75: \\(\\.\\([0-9]*\\)\\|\\)\
76: \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
77: "Regular expression to match floating point numbers. Extract matches:
78: 1 - minus sign
79: 2 - integer part
80: 4 - fractional part
81: 8 - minus sign for power of ten
82: 9 - power of ten
83: ")
84:
85: (defconst high-bit-mask (ash 1 maxbit)
86: "Masks all bits except the high-order (sign) bit.")
87:
88: (defconst second-bit-mask (ash 1 (1- maxbit))
89: "Masks all bits except the highest-order magnitude bit")
90:
91: ;; various useful floating point constants
92: (setq _f0 '(0 . 1))
93:
94: (setq _f1/2 '(4194304 . -23))
95:
96: (setq _f1 '(4194304 . -22))
97:
98: (setq _f10 '(5242880 . -19))
99:
100: ;; support for decimal conversion routines
101: (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
102: (aset powers-of-10 1 _f10)
103: (aset powers-of-10 2 '(6553600 . -16))
104: (aset powers-of-10 3 '(8192000 . -13))
105: (aset powers-of-10 4 '(5120000 . -9))
106: (aset powers-of-10 5 '(6400000 . -6))
107: (aset powers-of-10 6 '(8000000 . -3))
108:
109: (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
110: highest-power-of-10 (aref powers-of-10 decimal-digits))
111:
112: (defun fashl (fnum) ; floating-point arithmetic shift left
113: (cons (ash (car fnum) 1) (1- (cdr fnum))))
114:
115: (defun fashr (fnum) ; floating point arithmetic shift right
116: (cons (ash (car fnum) -1) (1+ (cdr fnum))))
117:
118: (defun normalize (fnum)
119: (if (> (car fnum) 0) ; make sure next-to-highest bit is set
120: (while (zerop (logand (car fnum) second-bit-mask))
121: (setq fnum (fashl fnum)))
122: (if (< (car fnum) 0) ; make sure highest bit is set
123: (while (zerop (logand (car fnum) high-bit-mask))
124: (setq fnum (fashl fnum)))
125: (setq fnum _f0))) ; "standard 0"
126: fnum)
127:
128: (defun abs (n) ; integer absolute value
129: (if (natnump n) n (- n)))
130:
131: (defun fabs (fnum) ; re-normalize after taking abs value
132: (normalize (cons (abs (car fnum)) (cdr fnum))))
133:
134: (defun xor (a b) ; logical exclusive or
135: (and (or a b) (not (and a b))))
136:
137: (defun same-sign (a b) ; two f-p numbers have same sign?
138: (not (xor (natnump (car a)) (natnump (car b)))))
139:
140: (defun extract-match (str i) ; used after string-match
141: (condition-case ()
142: (substring str (match-beginning i) (match-end i))
143: (error "")))
144:
145: ;; support for the multiplication function
146: (setq halfword-bits (/ mantissa-bits 2) ; bits in a halfword
147: masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
148: maskhi (lognot masklo) ; isolate the upper halfword
149: round-limit (ash 1 (/ halfword-bits 2)))
150:
151: (defun hihalf (n) ; return high halfword, shifted down
152: (ash (logand n maskhi) (- halfword-bits)))
153:
154: (defun lohalf (n) ; return low halfword
155: (logand n masklo))
156:
157: ;; Visible functions
158:
159: ;; Arithmetic functions
160: (defun f+ (a1 a2)
161: "Returns the sum of two floating point numbers."
162: (let ((f1 (fmax a1 a2))
163: (f2 (fmin a1 a2)))
164: (if (same-sign a1 a2)
165: (setq f1 (fashr f1) ; shift right to avoid overflow
166: f2 (fashr f2)))
167: (normalize
168: (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
169: (cdr f1)))))
170:
171: (defun f- (a1 &optional a2) ; unary or binary minus
172: "Returns the difference of two floating point numbers."
173: (if a2
174: (f+ a1 (f- a2))
175: (normalize (cons (- (car a1)) (cdr a1)))))
176:
177: (defun f* (a1 a2) ; multiply in halfword chunks
178: "Returns the product of two floating point numbers."
179: (let* ((i1 (car (fabs a1)))
180: (i2 (car (fabs a2)))
181: (sign (not (same-sign a1 a2)))
182: (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
183: (lohalf (* (hihalf i1) (lohalf i2)))
184: (lohalf (* (lohalf i1) (hihalf i2)))))
185: (prodhi (+ (* (hihalf i1) (hihalf i2))
186: (hihalf (* (hihalf i1) (lohalf i2)))
187: (hihalf (* (lohalf i1) (hihalf i2)))
188: (hihalf prodlo))))
189: (if (> (lohalf prodlo) round-limit)
190: (setq prodhi (1+ prodhi))) ; round off truncated bits
191: (normalize
192: (cons (if sign (- prodhi) prodhi)
193: (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
194:
195: (defun f/ (a1 a2) ; SLOW subtract-and-shift algorithm
196: "Returns the quotient of two floating point numbers."
197: (if (zerop (car a2)) ; if divide by 0
198: (signal 'arith-error (list "attempt to divide by zero" a1 a2))
199: (let ((bits (1- maxbit))
200: (quotient 0)
201: (dividend (car (fabs a1)))
202: (divisor (car (fabs a2)))
203: (sign (not (same-sign a1 a2))))
204: (while (natnump bits)
205: (if (< (- dividend divisor) 0)
206: (setq quotient (ash quotient 1))
207: (setq quotient (1+ (ash quotient 1))
208: dividend (- dividend divisor)))
209: (setq dividend (ash dividend 1)
210: bits (1- bits)))
211: (normalize
212: (cons (if sign (- quotient) quotient)
213: (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
214:
215: (defun f% (a1 a2)
216: "Returns the remainder of first floating point number divided by second."
217: (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
218:
219:
220: ;; Comparison functions
221: (defun f= (a1 a2)
222: "Returns t if two floating point numbers are equal, nil otherwise."
223: (equal a1 a2))
224:
225: (defun f> (a1 a2)
226: "Returns t if first floating point number is greater than second,
227: nil otherwise."
228: (cond ((and (natnump (car a1)) (< (car a2) 0))
229: t) ; a1 nonnegative, a2 negative
230: ((and (> (car a1) 0) (<= (car a2) 0))
231: t) ; a1 positive, a2 nonpositive
232: ((and (<= (car a1) 0) (natnump (car a2)))
233: nil) ; a1 nonpos, a2 nonneg
234: ((/= (cdr a1) (cdr a2)) ; same signs. exponents differ
235: (> (cdr a1) (cdr a2))) ; compare the mantissas.
236: (t
237: (> (car a1) (car a2))))) ; same exponents.
238:
239: (defun f>= (a1 a2)
240: "Returns t if first floating point number is greater than or equal to
241: second, nil otherwise."
242: (or (f> a1 a2) (f= a1 a2)))
243:
244: (defun f< (a1 a2)
245: "Returns t if first floating point number is less than second,
246: nil otherwise."
247: (not (f>= a1 a2)))
248:
249: (defun f<= (a1 a2)
250: "Returns t if first floating point number is less than or equal to
251: second, nil otherwise."
252: (not (f> a1 a2)))
253:
254: (defun f/= (a1 a2)
255: "Returns t if first floating point number is not equal to second,
256: nil otherwise."
257: (not (f= a1 a2)))
258:
259: (defun fmin (a1 a2)
260: "Returns the minimum of two floating point numbers."
261: (if (f< a1 a2) a1 a2))
262:
263: (defun fmax (a1 a2)
264: "Returns the maximum of two floating point numbers."
265: (if (f> a1 a2) a1 a2))
266:
267: (defun fzerop (fnum)
268: "Returns t if the floating point number is zero, nil otherwise."
269: (= (car fnum) 0))
270:
271: (defun floatp (fnum)
272: "Returns t if the arg is a floating point number, nil otherwise."
273: (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
274:
275: ;; Conversion routines
276: (defun f (int)
277: "Convert the integer argument to floating point, like a C cast operator."
278: (normalize (cons int '0)))
279:
280: (defun int-to-hex-string (int)
281: "Convert the integer argument to a C-style hexadecimal string."
282: (let ((shiftval -20)
283: (str "0x")
284: (hex-chars "0123456789ABCDEF"))
285: (while (<= shiftval 0)
286: (setq str (concat str (char-to-string
287: (aref hex-chars
288: (logand (lsh int shiftval) 15))))
289: shiftval (+ shiftval 4)))
290: str))
291:
292: (defun ftrunc (fnum) ; truncate fractional part
293: "Truncate the fractional part of a floating point number."
294: (cond ((natnump (cdr fnum)) ; it's all integer, return number as is
295: fnum)
296: ((<= (cdr fnum) (- maxbit)) ; it's all fractional, return 0
297: '(0 . 1))
298: (t ; otherwise mask out fractional bits
299: (let ((mant (car fnum)) (exp (cdr fnum)))
300: (normalize
301: (cons (if (natnump mant) ; if negative, use absolute value
302: (ash (ash mant exp) (- exp))
303: (- (ash (ash (- mant) exp) (- exp))))
304: exp))))))
305:
306: (defun fint (fnum) ; truncate and convert to integer
307: "Convert the floating point number to integer, with truncation,
308: like a C cast operator."
309: (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
310: (cond ((>= texp mantissa-bits) ; too high, return "maxint"
311: mantissa-maxval)
312: ((<= texp (- mantissa-bits)) ; too low, return "minint"
313: mantissa-minval)
314: (t ; in range
315: (ash tint texp))))) ; shift so that exponent is 0
316:
317: (defun float-to-string (fnum &optional sci)
318: "Convert the floating point number to a decimal string.
319: Optional second argument non-nil means use scientific notation."
320: (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
321: (power 0) (result 0) (str "")
322: (temp 0) (pow10 _f1))
323:
324: (if (f= fnum _f0)
325: "0"
326: (if (f>= value _f1) ; find largest power of 10 <= value
327: (progn ; value >= 1, power is positive
328: (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
329: (setq pow10 temp
330: power (+ power decimal-digits)))
331: (while (f<= (setq temp (f* pow10 _f10)) value)
332: (setq pow10 temp
333: power (1+ power))))
334: (progn ; value < 1, power is negative
335: (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
336: (setq pow10 temp
337: power (- power decimal-digits)))
338: (while (f> pow10 value)
339: (setq pow10 (f/ pow10 _f10)
340: power (1- power)))))
341: ; get value in range 100000 to 999999
342: (setq value (f* (f/ value pow10) all-decimal-digs-minval)
343: result (ftrunc value))
344: (if (f> (f- value result) _f1/2) ; round up if remainder > 0.5
345: (setq str (int-to-string (1+ (fint result))))
346: (setq str (int-to-string (fint result))))
347:
348: (if sci ; scientific notation
349: (setq str (concat (substring str 0 1) "." (substring str 1)
350: "E" (int-to-string power)))
351:
352: ; regular decimal string
353: (cond ((>= power (1- decimal-digits))
354: ; large power, append zeroes
355: (let ((zeroes (- power decimal-digits)))
356: (while (natnump zeroes)
357: (setq str (concat str "0")
358: zeroes (1- zeroes)))))
359:
360: ; negative power, prepend decimal
361: ((< power 0) ; point and zeroes
362: (let ((zeroes (- (- power) 2)))
363: (while (natnump zeroes)
364: (setq str (concat "0" str)
365: zeroes (1- zeroes)))
366: (setq str (concat "0." str))))
367:
368: (t ; in range, insert decimal point
369: (setq str (concat
370: (substring str 0 (1+ power))
371: "."
372: (substring str (1+ power)))))))
373:
374: (if sign ; if negative, prepend minus sign
375: (concat "-" str)
376: str))))
377:
378:
379: ;; string to float conversion.
380: ;; accepts scientific notation, but ignores anything after the first two
381: ;; digits of the exponent.
382: (defun string-to-float (str)
383: "Convert the string to a floating point number.
384: Accepts a decimal string in scientific notation,
385: with exponent preceded by either E or e.
386: Only the 6 most significant digits of the integer and fractional parts
387: are used; only the first two digits of the exponent are used.
388: Negative signs preceding both the decimal number and the exponent
389: are recognized."
390:
391: (if (string-match floating-point-regexp str 0)
392: (let (power)
393: (f*
394: ; calculate the mantissa
395: (let* ((int-subst (extract-match str 2))
396: (fract-subst (extract-match str 4))
397: (digit-string (concat int-subst fract-subst))
398: (mant-sign (equal (extract-match str 1) "-"))
399: (leading-0s 0) (round-up nil))
400:
401: ; get rid of leading 0's
402: (setq power (- (length int-subst) decimal-digits))
403: (while (and (< leading-0s (length digit-string))
404: (= (aref digit-string leading-0s) ?0))
405: (setq leading-0s (1+ leading-0s)))
406: (setq power (- power leading-0s)
407: digit-string (substring digit-string leading-0s))
408:
409: ; if more than 6 digits, round off
410: (if (> (length digit-string) decimal-digits)
411: (setq round-up (>= (aref digit-string decimal-digits) ?5)
412: digit-string (substring digit-string 0 decimal-digits))
413: (setq power (+ power (- decimal-digits (length digit-string)))))
414:
415: ; round up and add minus sign, if necessary
416: (f (* (+ (string-to-int digit-string)
417: (if round-up 1 0))
418: (if mant-sign -1 1))))
419:
420: ; calculate the exponent (power of ten)
421: (let* ((expt-subst (extract-match str 9))
422: (expt-sign (equal (extract-match str 8) "-"))
423: (expt 0) (chunks 0) (tens 0) (exponent _f1)
424: (func 'f*))
425:
426: (setq expt (+ (* (string-to-int
427: (substring expt-subst 0
428: (min expt-digits (length expt-subst))))
429: (if expt-sign -1 1))
430: power))
431: (if (< expt 0) ; if power of 10 negative
432: (setq expt (- expt) ; take abs val of exponent
433: func 'f/)) ; and set up to divide, not multiply
434:
435: (setq chunks (/ expt decimal-digits)
436: tens (% expt decimal-digits))
437: ; divide or multiply by "chunks" of 10**6
438: (while (> chunks 0)
439: (setq exponent (funcall func exponent highest-power-of-10)
440: chunks (1- chunks)))
441: ; divide or multiply by remaining power of ten
442: (funcall func exponent (aref powers-of-10 tens)))))
443:
444: _f0)) ; if invalid, return 0
445:
446:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.