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