summaryrefslogtreecommitdiff
path: root/lisp/emacs-lisp/float.el
blob: 85d9b4db78c358f44d400c76cfe5bfff947bc6e0 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
;;; float.el --- floating point arithmetic package.

;; Author: Bill Rosenblatt
;; Maintainer: FSF
;; Last-Modified: 16 Mar 1992

;; Copyright (C) 1986 Free Software Foundation, Inc.

;; This file is part of GNU Emacs.

;; GNU Emacs is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs; see the file COPYING.  If not, write to
;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.

;;; Commentary:

;; Floating point numbers are represented by dot-pairs (mant . exp)
;; where mant is the 24-bit signed integral mantissa and exp is the
;; base 2 exponent.
;;
;; Emacs LISP supports a 24-bit signed integer data type, which has a
;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
;; This gives six significant decimal digit accuracy.  Exponents can
;; be anything in the range -(2**23) to +(2**23)-1.
;;
;; User interface:
;; function f converts from integer to floating point
;; function string-to-float converts from string to floating point
;; function fint converts a floating point to integer (with truncation)
;; function float-to-string converts from floating point to string
;;                   
;; Caveats:
;; -  Exponents outside of the range of +/-100 or so will cause certain 
;;    functions (especially conversion routines) to take forever.
;; -  Very little checking is done for fixed point overflow/underflow.
;; -  No checking is done for over/underflow of the exponent
;;    (hardly necessary when exponent can be 2**23).
;; 
;;
;; Bill Rosenblatt
;; June 20, 1986
;;

;;; Code:

;; fundamental implementation constants
(defconst exp-base 2
  "Base of exponent in this floating point representation.")

(defconst mantissa-bits 24
  "Number of significant bits in this floating point representation.")

(defconst decimal-digits 6
  "Number of decimal digits expected to be accurate.")

(defconst expt-digits 2
  "Maximum permitted digits in a scientific notation exponent.")

;; other constants
(defconst maxbit (1- mantissa-bits)
  "Number of highest bit")

(defconst mantissa-maxval (1- (ash 1 maxbit))
  "Maximum permissable value of mantissa")

(defconst mantissa-minval (ash 1 maxbit)
  "Minimum permissable value of mantissa")

(defconst floating-point-regexp
  "^[ \t]*\\(-?\\)\\([0-9]*\\)\
\\(\\.\\([0-9]*\\)\\|\\)\
\\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  "Regular expression to match floating point numbers.  Extract matches:
1 - minus sign
2 - integer part
4 - fractional part
8 - minus sign for power of ten
9 - power of ten
")

(defconst high-bit-mask (ash 1 maxbit)
  "Masks all bits except the high-order (sign) bit.")

(defconst second-bit-mask (ash 1 (1- maxbit))
  "Masks all bits except the highest-order magnitude bit")

;; various useful floating point constants
(setq _f0 '(0 . 1))

(setq _f1/2 '(4194304 . -23))

(setq _f1 '(4194304 . -22))

(setq _f10 '(5242880 . -19))

;; support for decimal conversion routines
(setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
(aset powers-of-10 1 _f10)
(aset powers-of-10 2 '(6553600 . -16))
(aset powers-of-10 3 '(8192000 . -13))
(aset powers-of-10 4 '(5120000 . -9))
(aset powers-of-10 5 '(6400000 . -6))
(aset powers-of-10 6 '(8000000 . -3))

(setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
      highest-power-of-10 (aref powers-of-10 decimal-digits))

(defun fashl (fnum)			; floating-point arithmetic shift left
  (cons (ash (car fnum) 1) (1- (cdr fnum))))

(defun fashr (fnum)			; floating point arithmetic shift right
  (cons (ash (car fnum) -1) (1+ (cdr fnum))))

(defun normalize (fnum)
  (if (> (car fnum) 0)			; make sure next-to-highest bit is set
      (while (zerop (logand (car fnum) second-bit-mask))
	(setq fnum (fashl fnum)))
    (if (< (car fnum) 0)		; make sure highest bit is set
	(while (zerop (logand (car fnum) high-bit-mask))
	  (setq fnum (fashl fnum)))
      (setq fnum _f0)))			; "standard 0"
  fnum)
      
(defun abs (n)				; integer absolute value
  (if (>= n 0) n (- n)))

(defun fabs (fnum)			; re-normalize after taking abs value
  (normalize (cons (abs (car fnum)) (cdr fnum))))

(defun xor (a b)			; logical exclusive or
  (and (or a b) (not (and a b))))

(defun same-sign (a b)			; two f-p numbers have same sign?
  (not (xor (natnump (car a)) (natnump (car b)))))

(defun extract-match (str i)		; used after string-match
  (condition-case ()
      (substring str (match-beginning i) (match-end i))
    (error "")))

;; support for the multiplication function
(setq halfword-bits (/ mantissa-bits 2)	; bits in a halfword
      masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
      maskhi (lognot masklo)		; isolate the upper halfword
      round-limit (ash 1 (/ halfword-bits 2)))

(defun hihalf (n)			; return high halfword, shifted down
  (ash (logand n maskhi) (- halfword-bits)))

(defun lohalf (n)			; return low halfword
  (logand n masklo))

;; Visible functions

;; Arithmetic functions
(defun f+ (a1 a2)
  "Returns the sum of two floating point numbers."
  (let ((f1 (fmax a1 a2))
	(f2 (fmin a1 a2)))
    (if (same-sign a1 a2)
	(setq f1 (fashr f1)		; shift right to avoid overflow
	      f2 (fashr f2)))
    (normalize
     (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
	   (cdr f1)))))

(defun f- (a1 &optional a2)		; unary or binary minus
  "Returns the difference of two floating point numbers."
  (if a2
      (f+ a1 (f- a2))
    (normalize (cons (- (car a1)) (cdr a1)))))

(defun f* (a1 a2)			; multiply in halfword chunks
  "Returns the product of two floating point numbers."
  (let* ((i1 (car (fabs a1)))
	 (i2 (car (fabs a2)))
	 (sign (not (same-sign a1 a2)))
	 (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
		    (lohalf (* (hihalf i1) (lohalf i2)))
		    (lohalf (* (lohalf i1) (hihalf i2)))))
	 (prodhi (+ (* (hihalf i1) (hihalf i2))
		    (hihalf (* (hihalf i1) (lohalf i2)))
		    (hihalf (* (lohalf i1) (hihalf i2)))
		    (hihalf prodlo))))
    (if (> (lohalf prodlo) round-limit)
	(setq prodhi (1+ prodhi)))	; round off truncated bits
    (normalize
     (cons (if sign (- prodhi) prodhi)
	   (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))

(defun f/ (a1 a2)			; SLOW subtract-and-shift algorithm
  "Returns the quotient of two floating point numbers."
  (if (zerop (car a2))			; if divide by 0
      (signal 'arith-error (list "attempt to divide by zero" a1 a2))
    (let ((bits (1- maxbit))
	  (quotient 0) 
	  (dividend (car (fabs a1)))
	  (divisor (car (fabs a2)))
	  (sign (not (same-sign a1 a2))))
      (while (natnump bits)
	(if (< (- dividend divisor) 0)
	    (setq quotient (ash quotient 1))
	  (setq quotient (1+ (ash quotient 1))
		dividend (- dividend divisor)))
	(setq dividend (ash dividend 1)
	      bits (1- bits)))
      (normalize
       (cons (if sign (- quotient) quotient)
	     (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
  
(defun f% (a1 a2)
  "Returns the remainder of first floating point number divided by second."
  (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
	  

;; Comparison functions
(defun f= (a1 a2)
  "Returns t if two floating point numbers are equal, nil otherwise."
  (equal a1 a2))

(defun f> (a1 a2)
  "Returns t if first floating point number is greater than second,
nil otherwise."
  (cond ((and (natnump (car a1)) (< (car a2) 0)) 
	 t)				; a1 nonnegative, a2 negative
	((and (> (car a1) 0) (<= (car a2) 0))
	 t)				; a1 positive, a2 nonpositive
	((and (<= (car a1) 0) (natnump (car a2)))
	 nil)				; a1 nonpos, a2 nonneg
	((/= (cdr a1) (cdr a2))		; same signs.  exponents differ
	 (> (cdr a1) (cdr a2)))		; compare the mantissas.
	(t
	 (> (car a1) (car a2)))))	; same exponents.

(defun f>= (a1 a2)
  "Returns t if first floating point number is greater than or equal to 
second, nil otherwise."
  (or (f> a1 a2) (f= a1 a2)))

(defun f< (a1 a2)
  "Returns t if first floating point number is less than second,
nil otherwise."
  (not (f>= a1 a2)))

(defun f<= (a1 a2)
  "Returns t if first floating point number is less than or equal to
second, nil otherwise."
  (not (f> a1 a2)))

(defun f/= (a1 a2)
  "Returns t if first floating point number is not equal to second,
nil otherwise."
  (not (f= a1 a2)))

(defun fmin (a1 a2)
  "Returns the minimum of two floating point numbers."
  (if (f< a1 a2) a1 a2))

(defun fmax (a1 a2)
  "Returns the maximum of two floating point numbers."
  (if (f> a1 a2) a1 a2))
      
(defun fzerop (fnum)
  "Returns t if the floating point number is zero, nil otherwise."
  (= (car fnum) 0))

(defun floatp (fnum)
  "Returns t if the arg is a floating point number, nil otherwise."
  (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))

;; Conversion routines
(defun f (int)
  "Convert the integer argument to floating point, like a C cast operator."
  (normalize (cons int '0)))

(defun int-to-hex-string (int)
  "Convert the integer argument to a C-style hexadecimal string."
  (let ((shiftval -20)
	(str "0x")
	(hex-chars "0123456789ABCDEF"))
    (while (<= shiftval 0)
      (setq str (concat str (char-to-string 
			(aref hex-chars
			      (logand (lsh int shiftval) 15))))
	    shiftval (+ shiftval 4)))
    str))

(defun ftrunc (fnum)			; truncate fractional part
  "Truncate the fractional part of a floating point number."
  (cond ((natnump (cdr fnum))		; it's all integer, return number as is
	 fnum)
	((<= (cdr fnum) (- maxbit))	; it's all fractional, return 0
	 '(0 . 1))
	(t				; otherwise mask out fractional bits
	 (let ((mant (car fnum)) (exp (cdr fnum)))
	   (normalize 
	    (cons (if (natnump mant)	; if negative, use absolute value
		      (ash (ash mant exp) (- exp))
		    (- (ash (ash (- mant) exp) (- exp))))
		  exp))))))

(defun fint (fnum)			; truncate and convert to integer
  "Convert the floating point number to integer, with truncation, 
like a C cast operator."
  (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
    (cond ((>= texp mantissa-bits)	; too high, return "maxint"
	   mantissa-maxval)
	  ((<= texp (- mantissa-bits))	; too low, return "minint"
	   mantissa-minval)
	  (t				; in range
	   (ash tint texp)))))		; shift so that exponent is 0

(defun float-to-string (fnum &optional sci)
  "Convert the floating point number to a decimal string.
Optional second argument non-nil means use scientific notation."
  (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
	 (power 0) (result 0) (str "") 
	 (temp 0) (pow10 _f1))

    (if (f= fnum _f0)
	"0"
      (if (f>= value _f1)			; find largest power of 10 <= value
	  (progn				; value >= 1, power is positive
	    (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
	      (setq pow10 temp
		    power (+ power decimal-digits)))
	    (while (f<= (setq temp (f* pow10 _f10)) value)
	      (setq pow10 temp
		    power (1+ power))))
	(progn				; value < 1, power is negative
	  (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
	    (setq pow10 temp
		  power (- power decimal-digits)))
	  (while (f> pow10 value)
	    (setq pow10 (f/ pow10 _f10)
		  power (1- power)))))
					  ; get value in range 100000 to 999999
      (setq value (f* (f/ value pow10) all-decimal-digs-minval)
	    result (ftrunc value))
      (let (int)
	(if (f> (f- value result) _f1/2)	; round up if remainder > 0.5
	    (setq int (1+ (fint result)))
	  (setq int (fint result)))
	(setq str (int-to-string int))
	(if (>= int 1000000)
	    (setq power (1+ power))))

      (if sci				; scientific notation
	  (setq str (concat (substring str 0 1) "." (substring str 1)
			    "E" (int-to-string power)))

					  ; regular decimal string
	(cond ((>= power (1- decimal-digits))
					  ; large power, append zeroes
	       (let ((zeroes (- power decimal-digits)))
		 (while (natnump zeroes)
		   (setq str (concat str "0")
			 zeroes (1- zeroes)))))

					  ; negative power, prepend decimal
	      ((< power 0)		; point and zeroes
	       (let ((zeroes (- (- power) 2)))
		 (while (natnump zeroes)
		   (setq str (concat "0" str)
			 zeroes (1- zeroes)))
		 (setq str (concat "0." str))))

	      (t				; in range, insert decimal point
	       (setq str (concat
			  (substring str 0 (1+ power))
			  "."
			  (substring str (1+ power)))))))

      (if sign				; if negative, prepend minus sign
	  (concat "-" str)
	str))))

    
;; string to float conversion.
;; accepts scientific notation, but ignores anything after the first two
;; digits of the exponent.
(defun string-to-float (str)
  "Convert the string to a floating point number.
Accepts a decimal string in scientific notation,  with exponent preceded
by either E or e.  Only the six most significant digits of the integer
and fractional parts are used; only the first two digits of the exponent
are used.  Negative signs preceding both the decimal number and the exponent
are recognized."

  (if (string-match floating-point-regexp str 0)
      (let (power)
	(f*
	 ; calculate the mantissa
	 (let* ((int-subst (extract-match str 2))
		(fract-subst (extract-match str 4))
		(digit-string (concat int-subst fract-subst))
		(mant-sign (equal (extract-match str 1) "-"))
		(leading-0s 0) (round-up nil))

	   ; get rid of leading 0's
	   (setq power (- (length int-subst) decimal-digits))
	   (while (and (< leading-0s (length digit-string))
		       (= (aref digit-string leading-0s) ?0))
	     (setq leading-0s (1+ leading-0s)))
	   (setq power (- power leading-0s)
		 digit-string (substring digit-string leading-0s))
	   
	   ; if more than 6 digits, round off
	   (if (> (length digit-string) decimal-digits)
	       (setq round-up (>= (aref digit-string decimal-digits) ?5)
		     digit-string (substring digit-string 0 decimal-digits))
	     (setq power (+ power (- decimal-digits (length digit-string)))))

	   ; round up and add minus sign, if necessary
	   (f (* (+ (string-to-int digit-string)
		    (if round-up 1 0))
		 (if mant-sign -1 1))))
	   
	 ; calculate the exponent (power of ten)
	 (let* ((expt-subst (extract-match str 9))
		(expt-sign (equal (extract-match str 8) "-"))
		(expt 0) (chunks 0) (tens 0) (exponent _f1)
		(func 'f*))
 
	   (setq expt (+ (* (string-to-int
			     (substring expt-subst 0
					(min expt-digits (length expt-subst))))
			    (if expt-sign -1 1))
			 power))
	   (if (< expt 0)		; if power of 10 negative
	       (setq expt (- expt)	; take abs val of exponent
		     func 'f/))		; and set up to divide, not multiply

	   (setq chunks (/ expt decimal-digits)
		 tens (% expt decimal-digits))
	   ; divide or multiply by "chunks" of 10**6
	   (while (> chunks 0)	
	     (setq exponent (funcall func exponent highest-power-of-10)
		   chunks (1- chunks)))
	   ; divide or multiply by remaining power of ten
	   (funcall func exponent (aref powers-of-10 tens)))))
		  
    _f0))				; if invalid, return 0

(provide 'float)

;;; float.el ends here