1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/clojure/contrib/math.clj	Sat Aug 21 06:25:44 2010 -0400
     1.3 @@ -0,0 +1,247 @@
     1.4 +;;; math.clj: math functions that deal intelligently with the various
     1.5 +;;; types in Clojure's numeric tower, as well as math functions
     1.6 +;;; commonly found in Scheme implementations.
     1.7 +
     1.8 +;; by Mark Engelberg (mark.engelberg@gmail.com)
     1.9 +;; January 17, 2009
    1.10 +
    1.11 +;; expt - (expt x y) is x to the yth power, returns an exact number
    1.12 +;;   if the base is an exact number, and the power is an integer,
    1.13 +;;   otherwise returns a double.
    1.14 +;; abs - (abs n) is the absolute value of n
    1.15 +;; gcd - (gcd m n) returns the greatest common divisor of m and n
    1.16 +;; lcm - (lcm m n) returns the least common multiple of m and n
    1.17 +
    1.18 +;; The behavior of the next three functions on doubles is consistent
    1.19 +;; with the behavior of the corresponding functions
    1.20 +;; in Java's Math library, but on exact numbers, returns an integer.
    1.21 +
    1.22 +;; floor - (floor n) returns the greatest integer less than or equal to n.
    1.23 +;;   If n is an exact number, floor returns an integer,
    1.24 +;;   otherwise a double.
    1.25 +;; ceil - (ceil n) returns the least integer greater than or equal to n.
    1.26 +;;   If n is an exact number, ceil returns an integer,
    1.27 +;;   otherwise a double.
    1.28 +;; round - (round n) rounds to the nearest integer.
    1.29 +;;   round always returns an integer.  round rounds up for values
    1.30 +;;   exactly in between two integers.
    1.31 +
    1.32 +
    1.33 +;; sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,
    1.34 +;;   specifically, if the input is an exact number, and is a square
    1.35 +;;   of an exact number, the output will be exact.  The downside
    1.36 +;;   is that for the common case (inexact square root), some extra
    1.37 +;;   computation is done to look for an exact square root first.
    1.38 +;;   So if you need blazingly fast square root performance, and you
    1.39 +;;   know you're just going to need a double result, you're better
    1.40 +;;   off calling java's Math/sqrt, or alternatively, you could just
    1.41 +;;   convert your input to a double before calling this sqrt function.
    1.42 +;;   If Clojure ever gets complex numbers, then this function will
    1.43 +;;   need to be updated (so negative inputs yield complex outputs).
    1.44 +;; exact-integer-sqrt - Implements a math function from the R6RS Scheme
    1.45 +;;   standard.  (exact-integer-sqrt k) where k is a non-negative integer,
    1.46 +;;   returns [s r] where k = s^2+r and k < (s+1)^2.  In other words, it
    1.47 +;;   returns the floor of the square root and the "remainder".
    1.48 +
    1.49 +(ns 
    1.50 +  ^{:author "Mark Engelberg",
    1.51 +     :doc "Math functions that deal intelligently with the various
    1.52 +types in Clojure's numeric tower, as well as math functions
    1.53 +commonly found in Scheme implementations.
    1.54 +
    1.55 +expt - (expt x y) is x to the yth power, returns an exact number
    1.56 +  if the base is an exact number, and the power is an integer,
    1.57 +  otherwise returns a double.
    1.58 +abs - (abs n) is the absolute value of n
    1.59 +gcd - (gcd m n) returns the greatest common divisor of m and n
    1.60 +lcm - (lcm m n) returns the least common multiple of m and n
    1.61 +
    1.62 +The behavior of the next three functions on doubles is consistent
    1.63 +with the behavior of the corresponding functions
    1.64 +in Java's Math library, but on exact numbers, returns an integer.
    1.65 +
    1.66 +floor - (floor n) returns the greatest integer less than or equal to n.
    1.67 +  If n is an exact number, floor returns an integer,
    1.68 +  otherwise a double.
    1.69 +ceil - (ceil n) returns the least integer greater than or equal to n.
    1.70 +  If n is an exact number, ceil returns an integer,
    1.71 +  otherwise a double.
    1.72 +round - (round n) rounds to the nearest integer.
    1.73 +  round always returns an integer.  round rounds up for values
    1.74 +  exactly in between two integers.
    1.75 +
    1.76 +
    1.77 +sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,
    1.78 +  specifically, if the input is an exact number, and is a square
    1.79 +  of an exact number, the output will be exact.  The downside
    1.80 +  is that for the common case (inexact square root), some extra
    1.81 +  computation is done to look for an exact square root first.
    1.82 +  So if you need blazingly fast square root performance, and you
    1.83 +  know you're just going to need a double result, you're better
    1.84 +  off calling java's Math/sqrt, or alternatively, you could just
    1.85 +  convert your input to a double before calling this sqrt function.
    1.86 +  If Clojure ever gets complex numbers, then this function will
    1.87 +  need to be updated (so negative inputs yield complex outputs).
    1.88 +exact-integer-sqrt - Implements a math function from the R6RS Scheme
    1.89 +  standard.  (exact-integer-sqrt k) where k is a non-negative integer,
    1.90 +  returns [s r] where k = s^2+r and k < (s+1)^2.  In other words, it
    1.91 +  returns the floor of the square root and the "remainder".
    1.92 +"}
    1.93 +  clojure.contrib.math)
    1.94 +
    1.95 +(derive ::integer ::exact)
    1.96 +(derive java.lang.Integer ::integer)
    1.97 +(derive java.math.BigInteger ::integer)
    1.98 +(derive java.lang.Long ::integer)
    1.99 +(derive java.math.BigDecimal ::exact)
   1.100 +(derive clojure.lang.Ratio ::exact)
   1.101 +(derive java.lang.Double ::inexact)
   1.102 +(derive java.lang.Float ::inexact)
   1.103 +
   1.104 +(defmulti ^{:arglists '([base pow])
   1.105 +	     :doc "(expt base pow) is base to the pow power.
   1.106 +Returns an exact number if the base is an exact number and the power is an integer, otherwise returns a double."}
   1.107 +  expt (fn [x y] [(class x) (class y)]))
   1.108 +
   1.109 +(defn- expt-int [base pow]
   1.110 +  (loop [n pow, y (num 1), z base]
   1.111 +    (let [t (bit-and n 1), n (bit-shift-right n 1)]
   1.112 +      (cond
   1.113 +       (zero? t) (recur n y (* z z))
   1.114 +       (zero? n) (* z y)
   1.115 +       :else (recur n (* z y) (* z z))))))
   1.116 +
   1.117 +(defmethod expt [::exact ::integer] [base pow]
   1.118 +  (cond
   1.119 +   (pos? pow) (expt-int base pow)
   1.120 +   (zero? pow) 1
   1.121 +   :else (/ 1 (expt-int base (- pow)))))
   1.122 +
   1.123 +(defmethod expt :default [base pow] (Math/pow base pow))
   1.124 +
   1.125 +(defn abs "(abs n) is the absolute value of n" [n]
   1.126 +  (cond
   1.127 +   (not (number? n)) (throw (IllegalArgumentException.
   1.128 +			     "abs requires a number"))
   1.129 +   (neg? n) (- n)
   1.130 +   :else n))
   1.131 +
   1.132 +(defmulti ^{:arglists '([n])
   1.133 +	     :doc "(floor n) returns the greatest integer less than or equal to n.
   1.134 +If n is an exact number, floor returns an integer, otherwise a double."}
   1.135 +  floor class)
   1.136 +(defmethod floor ::integer [n] n)
   1.137 +(defmethod floor java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_FLOOR) (toBigInteger)))
   1.138 +(defmethod floor clojure.lang.Ratio [n]
   1.139 +  (if (pos? n) (quot (. n numerator) (. n denominator))
   1.140 +      (dec (quot (. n numerator) (. n denominator)))))
   1.141 +(defmethod floor :default [n]
   1.142 +  (Math/floor n))
   1.143 +
   1.144 +(defmulti ^{:arglists '([n])
   1.145 +	     :doc "(ceil n) returns the least integer greater than or equal to n.
   1.146 +If n is an exact number, ceil returns an integer, otherwise a double."}
   1.147 +  ceil class)
   1.148 +(defmethod ceil ::integer [n] n)
   1.149 +(defmethod ceil java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_CEILING) (toBigInteger)))
   1.150 +(defmethod ceil clojure.lang.Ratio [n]
   1.151 +  (if (pos? n) (inc (quot (. n numerator) (. n denominator)))
   1.152 +      (quot (. n numerator) (. n denominator))))
   1.153 +(defmethod ceil :default [n]
   1.154 +  (Math/ceil n))
   1.155 +
   1.156 +(defmulti ^{:arglists '([n])
   1.157 +	     :doc "(round n) rounds to the nearest integer.
   1.158 +round always returns an integer.  Rounds up for values exactly in between two integers."}
   1.159 +  round class)
   1.160 +(defmethod round ::integer [n] n)
   1.161 +(defmethod round java.math.BigDecimal [n] (floor (+ n 0.5M)))
   1.162 +(defmethod round clojure.lang.Ratio [n] (floor (+ n 1/2)))
   1.163 +(defmethod round :default [n] (Math/round n))
   1.164 +
   1.165 +(defn gcd "(gcd a b) returns the greatest common divisor of a and b" [a b]
   1.166 +  (if (or (not (integer? a)) (not (integer? b)))
   1.167 +    (throw (IllegalArgumentException. "gcd requires two integers"))  
   1.168 +    (loop [a (abs a) b (abs b)]
   1.169 +      (if (zero? b) a,
   1.170 +	  (recur b (mod a b))))))
   1.171 +
   1.172 +(defn lcm
   1.173 +  "(lcm a b) returns the least common multiple of a and b"
   1.174 +  [a b]
   1.175 +  (when (or (not (integer? a)) (not (integer? b)))
   1.176 +    (throw (IllegalArgumentException. "lcm requires two integers")))
   1.177 +  (cond (zero? a) 0
   1.178 +        (zero? b) 0
   1.179 +        :else (abs (* b (quot a (gcd a b))))))
   1.180 +
   1.181 +; Length of integer in binary, used as helper function for sqrt.
   1.182 +(defmulti ^{:private true} integer-length class)
   1.183 +(defmethod integer-length java.lang.Integer [n]
   1.184 +  (count (Integer/toBinaryString n)))
   1.185 +(defmethod integer-length java.lang.Long [n]
   1.186 +  (count (Long/toBinaryString n)))
   1.187 +(defmethod integer-length java.math.BigInteger [n]
   1.188 +  (count (. n toString 2)))
   1.189 +
   1.190 +;; Produces the largest integer less than or equal to the square root of n
   1.191 +;; Input n must be a non-negative integer
   1.192 +(defn- integer-sqrt [n]
   1.193 +  (cond
   1.194 +   (> n 24)
   1.195 +   (let [n-len (integer-length n)]
   1.196 +     (loop [init-value (if (even? n-len)
   1.197 +			 (bit-shift-left 1 (bit-shift-right n-len 1))
   1.198 +			 (bit-shift-left 2 (bit-shift-right n-len 1)))]
   1.199 +       (let [iterated-value (bit-shift-right (+ init-value (quot n init-value)) 1)]
   1.200 +	 (if (>= iterated-value init-value)
   1.201 +	   init-value
   1.202 +	   (recur iterated-value)))))
   1.203 +   (> n 15) 4
   1.204 +   (> n  8) 3
   1.205 +   (> n  3) 2
   1.206 +   (> n  0) 1
   1.207 +   (> n -1) 0))
   1.208 +
   1.209 +(defn exact-integer-sqrt "(exact-integer-sqrt n) expects a non-negative integer n, and returns [s r] where n = s^2+r and n < (s+1)^2.  In other words, it returns the floor of the square root and the 'remainder'.
   1.210 +For example, (exact-integer-sqrt 15) is [3 6] because 15 = 3^2+6."
   1.211 +  [n]
   1.212 +  (if (or (not (integer? n)) (neg? n))
   1.213 +    (throw (IllegalArgumentException. "exact-integer-sqrt requires a non-negative integer"))
   1.214 +    (let [isqrt (integer-sqrt n),
   1.215 +	  error (- n (* isqrt isqrt))]
   1.216 +      [isqrt error])))
   1.217 +
   1.218 +(defmulti ^{:arglists '([n])
   1.219 +	     :doc "Square root, but returns exact number if possible."}
   1.220 +  sqrt class)
   1.221 +(defmethod sqrt ::integer [n]
   1.222 +  (if (neg? n) Double/NaN
   1.223 +      (let [isqrt (integer-sqrt n),
   1.224 +	    error (- n (* isqrt isqrt))]
   1.225 +	(if (zero? error) isqrt
   1.226 +	    (Math/sqrt n)))))
   1.227 +
   1.228 +(defmethod sqrt clojure.lang.Ratio [n]
   1.229 +  (if (neg? n) Double/NaN
   1.230 +      (let [numerator (.numerator n),
   1.231 +	    denominator (.denominator n),
   1.232 +	    sqrtnum (sqrt numerator)]
   1.233 +	(if (float? sqrtnum)
   1.234 +	  (Math/sqrt n)
   1.235 +	  (let [sqrtden (sqrt denominator)]
   1.236 +	    (if (float? sqrtnum)
   1.237 +	      (Math/sqrt n)
   1.238 +	      (/ sqrtnum sqrtden)))))))
   1.239 +
   1.240 +(defmethod sqrt java.math.BigDecimal [n]
   1.241 +  (if (neg? n) Double/NaN
   1.242 +      (let [frac (rationalize n),
   1.243 +	    sqrtfrac (sqrt frac)]
   1.244 +	(if (ratio? sqrtfrac)
   1.245 +	  (/ (BigDecimal. (.numerator sqrtfrac))
   1.246 +	     (BigDecimal. (.denominator sqrtfrac)))
   1.247 +	  sqrtfrac))))
   1.248 +
   1.249 +(defmethod sqrt :default [n]
   1.250 +  (Math/sqrt n))