home *** CD-ROM | disk | FTP | other *** search
/ OpenStep 4.2J (Developer) / os42jdev.iso / NextDeveloper / Source / GNU / emacs / lisp / float.el < prev    next >
Lisp/Scheme  |  1992-02-01  |  15KB  |  460 lines
  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.  
  460.