(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.")
(defconst maxbit (1- mantissa-bits)
"Number of highest bit")
(defconst mantissa-maxval (1- (ash 1 maxbit))
"Maximum permissible value of mantissa")
(defconst mantissa-minval (ash 1 maxbit)
"Minimum permissible 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")
(defconst _f0 '(0 . 1))
(defconst _f1/2 '(4194304 . -23))
(defconst _f1 '(4194304 . -22))
(defconst _f10 '(5242880 . -19))
(defvar 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))
(defconst all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits)))
(defconst highest-power-of-10 (aref powers-of-10 decimal-digits))
(defun fashl (fnum) (cons (ash (car fnum) 1) (1- (cdr fnum))))
(defun fashr (fnum) (cons (ash (car fnum) -1) (1+ (cdr fnum))))
(defun normalize (fnum)
(if (> (car fnum) 0) (while (zerop (logand (car fnum) second-bit-mask))
(setq fnum (fashl fnum)))
(if (< (car fnum) 0) (while (zerop (logand (car fnum) high-bit-mask))
(setq fnum (fashl fnum)))
(setq fnum _f0))) fnum)
(defun abs (n) (if (>= n 0) n (- n)))
(defun fabs (fnum) (normalize (cons (abs (car fnum)) (cdr fnum))))
(defun xor (a b) (and (or a b) (not (and a b))))
(defun same-sign (a b) (not (xor (natnump (car a)) (natnump (car b)))))
(defun extract-match (str i) (condition-case ()
(substring str (match-beginning i) (match-end i))
(error "")))
(defconst halfword-bits (/ mantissa-bits 2)) (defconst masklo (1- (ash 1 halfword-bits))) (defconst maskhi (lognot masklo)) (defconst round-limit (ash 1 (/ halfword-bits 2)))
(defun hihalf (n) (ash (logand n maskhi) (- halfword-bits)))
(defun lohalf (n) (logand n masklo))
(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) f2 (fashr f2)))
(normalize
(cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
(cdr f1)))))
(defun f- (a1 &optional a2) "Returns the difference of two floating point numbers."
(if a2
(f+ a1 (f- a2))
(normalize (cons (- (car a1)) (cdr a1)))))
(defun f* (a1 a2) "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))) (normalize
(cons (if sign (- prodhi) prodhi)
(+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
(defun f/ (a1 a2) "Returns the quotient of two floating point numbers."
(if (zerop (car a2)) (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)))
(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) ((and (> (car a1) 0) (<= (car a2) 0))
t) ((and (<= (car a1) 0) (natnump (car a2)))
nil) ((/= (cdr a1) (cdr a2)) (> (cdr a1) (cdr a2))) (t
(> (car a1) (car a2)))))
(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))))
(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 the fractional part of a floating point number."
(cond ((natnump (cdr fnum)) fnum)
((<= (cdr fnum) (- maxbit)) '(0 . 1))
(t (let ((mant (car fnum)) (exp (cdr fnum)))
(normalize
(cons (if (natnump mant) (ash (ash mant exp) (- exp))
(- (ash (ash (- mant) exp) (- exp))))
exp))))))
(defun fint (fnum) "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) mantissa-maxval)
((<= texp (- mantissa-bits)) mantissa-minval)
(t (ash tint texp)))))
(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) (progn (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 (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)))))
(setq value (f* (f/ value pow10) all-decimal-digs-minval)
result (ftrunc value))
(let (int)
(if (f> (f- value result) _f1/2) (setq int (1+ (fint result)))
(setq int (fint result)))
(setq str (int-to-string int))
(if (>= int 1000000)
(setq power (1+ power))))
(if sci (setq str (concat (substring str 0 1) "." (substring str 1)
"E" (int-to-string power)))
(cond ((>= power (1- decimal-digits))
(let ((zeroes (- power decimal-digits)))
(while (natnump zeroes)
(setq str (concat str "0")
zeroes (1- zeroes)))))
((< power 0) (let ((zeroes (- (- power) 2)))
(while (natnump zeroes)
(setq str (concat "0" str)
zeroes (1- zeroes)))
(setq str (concat "0." str))))
(t (setq str (concat
(substring str 0 (1+ power))
"."
(substring str (1+ power)))))))
(if sign (concat "-" str)
str))))
(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*
(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))
(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 (> (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)))))
(f (* (+ (string-to-number digit-string)
(if round-up 1 0))
(if mant-sign -1 1))))
(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-number
(substring expt-subst 0
(min expt-digits (length expt-subst))))
(if expt-sign -1 1))
power))
(if (< expt 0) (setq expt (- expt) func 'f/))
(setq chunks (/ expt decimal-digits)
tens (% expt decimal-digits))
(while (> chunks 0)
(setq exponent (funcall func exponent highest-power-of-10)
chunks (1- chunks)))
(funcall func exponent (aref powers-of-10 tens)))))
_f0))
(provide 'float)