rlm@10: ;;; math.clj: math functions that deal intelligently with the various
rlm@10: ;;; types in Clojure's numeric tower, as well as math functions
rlm@10: ;;; commonly found in Scheme implementations.
rlm@10: 
rlm@10: ;; by Mark Engelberg (mark.engelberg@gmail.com)
rlm@10: ;; January 17, 2009
rlm@10: 
rlm@10: ;; expt - (expt x y) is x to the yth power, returns an exact number
rlm@10: ;;   if the base is an exact number, and the power is an integer,
rlm@10: ;;   otherwise returns a double.
rlm@10: ;; abs - (abs n) is the absolute value of n
rlm@10: ;; gcd - (gcd m n) returns the greatest common divisor of m and n
rlm@10: ;; lcm - (lcm m n) returns the least common multiple of m and n
rlm@10: 
rlm@10: ;; The behavior of the next three functions on doubles is consistent
rlm@10: ;; with the behavior of the corresponding functions
rlm@10: ;; in Java's Math library, but on exact numbers, returns an integer.
rlm@10: 
rlm@10: ;; floor - (floor n) returns the greatest integer less than or equal to n.
rlm@10: ;;   If n is an exact number, floor returns an integer,
rlm@10: ;;   otherwise a double.
rlm@10: ;; ceil - (ceil n) returns the least integer greater than or equal to n.
rlm@10: ;;   If n is an exact number, ceil returns an integer,
rlm@10: ;;   otherwise a double.
rlm@10: ;; round - (round n) rounds to the nearest integer.
rlm@10: ;;   round always returns an integer.  round rounds up for values
rlm@10: ;;   exactly in between two integers.
rlm@10: 
rlm@10: 
rlm@10: ;; sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,
rlm@10: ;;   specifically, if the input is an exact number, and is a square
rlm@10: ;;   of an exact number, the output will be exact.  The downside
rlm@10: ;;   is that for the common case (inexact square root), some extra
rlm@10: ;;   computation is done to look for an exact square root first.
rlm@10: ;;   So if you need blazingly fast square root performance, and you
rlm@10: ;;   know you're just going to need a double result, you're better
rlm@10: ;;   off calling java's Math/sqrt, or alternatively, you could just
rlm@10: ;;   convert your input to a double before calling this sqrt function.
rlm@10: ;;   If Clojure ever gets complex numbers, then this function will
rlm@10: ;;   need to be updated (so negative inputs yield complex outputs).
rlm@10: ;; exact-integer-sqrt - Implements a math function from the R6RS Scheme
rlm@10: ;;   standard.  (exact-integer-sqrt k) where k is a non-negative integer,
rlm@10: ;;   returns [s r] where k = s^2+r and k < (s+1)^2.  In other words, it
rlm@10: ;;   returns the floor of the square root and the "remainder".
rlm@10: 
rlm@10: (ns 
rlm@10:   ^{:author "Mark Engelberg",
rlm@10:      :doc "Math functions that deal intelligently with the various
rlm@10: types in Clojure's numeric tower, as well as math functions
rlm@10: commonly found in Scheme implementations.
rlm@10: 
rlm@10: expt - (expt x y) is x to the yth power, returns an exact number
rlm@10:   if the base is an exact number, and the power is an integer,
rlm@10:   otherwise returns a double.
rlm@10: abs - (abs n) is the absolute value of n
rlm@10: gcd - (gcd m n) returns the greatest common divisor of m and n
rlm@10: lcm - (lcm m n) returns the least common multiple of m and n
rlm@10: 
rlm@10: The behavior of the next three functions on doubles is consistent
rlm@10: with the behavior of the corresponding functions
rlm@10: in Java's Math library, but on exact numbers, returns an integer.
rlm@10: 
rlm@10: floor - (floor n) returns the greatest integer less than or equal to n.
rlm@10:   If n is an exact number, floor returns an integer,
rlm@10:   otherwise a double.
rlm@10: ceil - (ceil n) returns the least integer greater than or equal to n.
rlm@10:   If n is an exact number, ceil returns an integer,
rlm@10:   otherwise a double.
rlm@10: round - (round n) rounds to the nearest integer.
rlm@10:   round always returns an integer.  round rounds up for values
rlm@10:   exactly in between two integers.
rlm@10: 
rlm@10: 
rlm@10: sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,
rlm@10:   specifically, if the input is an exact number, and is a square
rlm@10:   of an exact number, the output will be exact.  The downside
rlm@10:   is that for the common case (inexact square root), some extra
rlm@10:   computation is done to look for an exact square root first.
rlm@10:   So if you need blazingly fast square root performance, and you
rlm@10:   know you're just going to need a double result, you're better
rlm@10:   off calling java's Math/sqrt, or alternatively, you could just
rlm@10:   convert your input to a double before calling this sqrt function.
rlm@10:   If Clojure ever gets complex numbers, then this function will
rlm@10:   need to be updated (so negative inputs yield complex outputs).
rlm@10: exact-integer-sqrt - Implements a math function from the R6RS Scheme
rlm@10:   standard.  (exact-integer-sqrt k) where k is a non-negative integer,
rlm@10:   returns [s r] where k = s^2+r and k < (s+1)^2.  In other words, it
rlm@10:   returns the floor of the square root and the "remainder".
rlm@10: "}
rlm@10:   clojure.contrib.math)
rlm@10: 
rlm@10: (derive ::integer ::exact)
rlm@10: (derive java.lang.Integer ::integer)
rlm@10: (derive java.math.BigInteger ::integer)
rlm@10: (derive java.lang.Long ::integer)
rlm@10: (derive java.math.BigDecimal ::exact)
rlm@10: (derive clojure.lang.Ratio ::exact)
rlm@10: (derive java.lang.Double ::inexact)
rlm@10: (derive java.lang.Float ::inexact)
rlm@10: 
rlm@10: (defmulti ^{:arglists '([base pow])
rlm@10: 	     :doc "(expt base pow) is base to the pow power.
rlm@10: Returns an exact number if the base is an exact number and the power is an integer, otherwise returns a double."}
rlm@10:   expt (fn [x y] [(class x) (class y)]))
rlm@10: 
rlm@10: (defn- expt-int [base pow]
rlm@10:   (loop [n pow, y (num 1), z base]
rlm@10:     (let [t (bit-and n 1), n (bit-shift-right n 1)]
rlm@10:       (cond
rlm@10:        (zero? t) (recur n y (* z z))
rlm@10:        (zero? n) (* z y)
rlm@10:        :else (recur n (* z y) (* z z))))))
rlm@10: 
rlm@10: (defmethod expt [::exact ::integer] [base pow]
rlm@10:   (cond
rlm@10:    (pos? pow) (expt-int base pow)
rlm@10:    (zero? pow) 1
rlm@10:    :else (/ 1 (expt-int base (- pow)))))
rlm@10: 
rlm@10: (defmethod expt :default [base pow] (Math/pow base pow))
rlm@10: 
rlm@10: (defn abs "(abs n) is the absolute value of n" [n]
rlm@10:   (cond
rlm@10:    (not (number? n)) (throw (IllegalArgumentException.
rlm@10: 			     "abs requires a number"))
rlm@10:    (neg? n) (- n)
rlm@10:    :else n))
rlm@10: 
rlm@10: (defmulti ^{:arglists '([n])
rlm@10: 	     :doc "(floor n) returns the greatest integer less than or equal to n.
rlm@10: If n is an exact number, floor returns an integer, otherwise a double."}
rlm@10:   floor class)
rlm@10: (defmethod floor ::integer [n] n)
rlm@10: (defmethod floor java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_FLOOR) (toBigInteger)))
rlm@10: (defmethod floor clojure.lang.Ratio [n]
rlm@10:   (if (pos? n) (quot (. n numerator) (. n denominator))
rlm@10:       (dec (quot (. n numerator) (. n denominator)))))
rlm@10: (defmethod floor :default [n]
rlm@10:   (Math/floor n))
rlm@10: 
rlm@10: (defmulti ^{:arglists '([n])
rlm@10: 	     :doc "(ceil n) returns the least integer greater than or equal to n.
rlm@10: If n is an exact number, ceil returns an integer, otherwise a double."}
rlm@10:   ceil class)
rlm@10: (defmethod ceil ::integer [n] n)
rlm@10: (defmethod ceil java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_CEILING) (toBigInteger)))
rlm@10: (defmethod ceil clojure.lang.Ratio [n]
rlm@10:   (if (pos? n) (inc (quot (. n numerator) (. n denominator)))
rlm@10:       (quot (. n numerator) (. n denominator))))
rlm@10: (defmethod ceil :default [n]
rlm@10:   (Math/ceil n))
rlm@10: 
rlm@10: (defmulti ^{:arglists '([n])
rlm@10: 	     :doc "(round n) rounds to the nearest integer.
rlm@10: round always returns an integer.  Rounds up for values exactly in between two integers."}
rlm@10:   round class)
rlm@10: (defmethod round ::integer [n] n)
rlm@10: (defmethod round java.math.BigDecimal [n] (floor (+ n 0.5M)))
rlm@10: (defmethod round clojure.lang.Ratio [n] (floor (+ n 1/2)))
rlm@10: (defmethod round :default [n] (Math/round n))
rlm@10: 
rlm@10: (defn gcd "(gcd a b) returns the greatest common divisor of a and b" [a b]
rlm@10:   (if (or (not (integer? a)) (not (integer? b)))
rlm@10:     (throw (IllegalArgumentException. "gcd requires two integers"))  
rlm@10:     (loop [a (abs a) b (abs b)]
rlm@10:       (if (zero? b) a,
rlm@10: 	  (recur b (mod a b))))))
rlm@10: 
rlm@10: (defn lcm
rlm@10:   "(lcm a b) returns the least common multiple of a and b"
rlm@10:   [a b]
rlm@10:   (when (or (not (integer? a)) (not (integer? b)))
rlm@10:     (throw (IllegalArgumentException. "lcm requires two integers")))
rlm@10:   (cond (zero? a) 0
rlm@10:         (zero? b) 0
rlm@10:         :else (abs (* b (quot a (gcd a b))))))
rlm@10: 
rlm@10: ; Length of integer in binary, used as helper function for sqrt.
rlm@10: (defmulti ^{:private true} integer-length class)
rlm@10: (defmethod integer-length java.lang.Integer [n]
rlm@10:   (count (Integer/toBinaryString n)))
rlm@10: (defmethod integer-length java.lang.Long [n]
rlm@10:   (count (Long/toBinaryString n)))
rlm@10: (defmethod integer-length java.math.BigInteger [n]
rlm@10:   (count (. n toString 2)))
rlm@10: 
rlm@10: ;; Produces the largest integer less than or equal to the square root of n
rlm@10: ;; Input n must be a non-negative integer
rlm@10: (defn- integer-sqrt [n]
rlm@10:   (cond
rlm@10:    (> n 24)
rlm@10:    (let [n-len (integer-length n)]
rlm@10:      (loop [init-value (if (even? n-len)
rlm@10: 			 (bit-shift-left 1 (bit-shift-right n-len 1))
rlm@10: 			 (bit-shift-left 2 (bit-shift-right n-len 1)))]
rlm@10:        (let [iterated-value (bit-shift-right (+ init-value (quot n init-value)) 1)]
rlm@10: 	 (if (>= iterated-value init-value)
rlm@10: 	   init-value
rlm@10: 	   (recur iterated-value)))))
rlm@10:    (> n 15) 4
rlm@10:    (> n  8) 3
rlm@10:    (> n  3) 2
rlm@10:    (> n  0) 1
rlm@10:    (> n -1) 0))
rlm@10: 
rlm@10: (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'.
rlm@10: For example, (exact-integer-sqrt 15) is [3 6] because 15 = 3^2+6."
rlm@10:   [n]
rlm@10:   (if (or (not (integer? n)) (neg? n))
rlm@10:     (throw (IllegalArgumentException. "exact-integer-sqrt requires a non-negative integer"))
rlm@10:     (let [isqrt (integer-sqrt n),
rlm@10: 	  error (- n (* isqrt isqrt))]
rlm@10:       [isqrt error])))
rlm@10: 
rlm@10: (defmulti ^{:arglists '([n])
rlm@10: 	     :doc "Square root, but returns exact number if possible."}
rlm@10:   sqrt class)
rlm@10: (defmethod sqrt ::integer [n]
rlm@10:   (if (neg? n) Double/NaN
rlm@10:       (let [isqrt (integer-sqrt n),
rlm@10: 	    error (- n (* isqrt isqrt))]
rlm@10: 	(if (zero? error) isqrt
rlm@10: 	    (Math/sqrt n)))))
rlm@10: 
rlm@10: (defmethod sqrt clojure.lang.Ratio [n]
rlm@10:   (if (neg? n) Double/NaN
rlm@10:       (let [numerator (.numerator n),
rlm@10: 	    denominator (.denominator n),
rlm@10: 	    sqrtnum (sqrt numerator)]
rlm@10: 	(if (float? sqrtnum)
rlm@10: 	  (Math/sqrt n)
rlm@10: 	  (let [sqrtden (sqrt denominator)]
rlm@10: 	    (if (float? sqrtnum)
rlm@10: 	      (Math/sqrt n)
rlm@10: 	      (/ sqrtnum sqrtden)))))))
rlm@10: 
rlm@10: (defmethod sqrt java.math.BigDecimal [n]
rlm@10:   (if (neg? n) Double/NaN
rlm@10:       (let [frac (rationalize n),
rlm@10: 	    sqrtfrac (sqrt frac)]
rlm@10: 	(if (ratio? sqrtfrac)
rlm@10: 	  (/ (BigDecimal. (.numerator sqrtfrac))
rlm@10: 	     (BigDecimal. (.denominator sqrtfrac)))
rlm@10: 	  sqrtfrac))))
rlm@10: 
rlm@10: (defmethod sqrt :default [n]
rlm@10:   (Math/sqrt n))