Added math.clj, which contains several math functions that are considered standard in Scheme implementations, adapted to Clojure's numeric tower when relevant, using Java's Math library for doubles.

Also added a math folder, which contains tests.clj (test cases for math.clj).

2 files changed, 284 insertions, 0 deletions

diff --git a/src/clojure/contrib/math.clj b/src/clojure/contrib/math.clj
new file mode 100644
index 00000000..787deeb2
--- /dev/null
+++ b/src/clojure/contrib/math.clj
@@ -0,0 +1,188 @@
+;;; math.clj: math functions that deal intelligently with the various

+;;; types in Clojure's numeric tower, as well as math functions

+;;; commonly found in Scheme implementations.

+

+;; by Mark Engelberg (mark.engelberg@gmail.com)

+;; January 17, 2009

+

+;; expt - (expt x y) is x to the yth power, returns an exact number

+;;   if the base is an exact number, and the power is an integer,

+;;   otherwise returns a double.

+;; abs - (abs n) is the absolute value of n

+;; gcd - (gcd m n) returns the greatest common divisor of m and n

+

+;; The behavior of the next three functions on doubles is consistent

+;; with the behavior of the corresponding functions

+;; in Java's Math library, but on exact numbers, returns an integer.

+

+;; floor - (floor n) returns the greatest integer less than or equal to n.

+;;   If n is an exact number, floor returns an integer,

+;;   otherwise a double.

+;; ceil - (ceil n) returns the least integer greater than or equal to n.

+;;   If n is an exact number, ceil returns an integer,

+;;   otherwise a double.

+;; round - (round n) rounds to the nearest integer.

+;;   round always returns an integer.  round rounds up for values

+;;   exactly in between two integers.

+

+

+;; sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,

+;;   specifically, if the input is an exact number, and is a square

+;;   of an exact number, the output will be exact.  The downside

+;;   is that for the common case (inexact square root), some extra

+;;   computation is done to look for an exact square root first.

+;;   So if you need blazingly fast square root performance, and you

+;;   know you're just going to need a double result, you're better

+;;   off calling java's Math/sqrt, or alternatively, you could just

+;;   convert your input to a double before calling this sqrt function.

+;;   If Clojure ever gets complex numbers, then this function will

+;;   need to be updated (so negative inputs yield complex outputs).

+;; exact-integer-sqrt - Implements a math function from the R6RS Scheme

+;;   standard.  (exact-integer-sqrt k) where k is a non-negative integer,

+;;   returns [s r] where k = s^2+r and k < (s+1)^2.  In other words, it

+;;   returns the floor of the square root and the "remainder".

+

+(ns clojure.contrib.math)

+

+(derive ::integer ::exact)

+(derive java.lang.Integer ::integer)

+(derive java.math.BigInteger ::integer)

+(derive java.lang.Long ::integer)

+(derive java.math.BigDecimal ::exact)

+(derive clojure.lang.Ratio ::exact)

+(derive java.lang.Double ::inexact)

+(derive java.lang.Float ::inexact)

+

+(defmulti expt (fn [x y] [(class x) (class y)]))

+

+(defn- expt-int [base pow]

+  (loop [n pow, y 1, z base]

+    (let [t (bit-and n 1), n (bit-shift-right n 1)]

+      (cond

+       (zero? t) (recur n y (* z z))

+       (zero? n) (* z y)

+       :else (recur n (* z y) (* z z))))))

+

+(defmethod expt [::exact ::integer] [base pow]

+  (cond

+   (pos? pow) (expt-int base pow)

+   (zero? pow) 1

+   :else (/ 1 (expt-int base (- pow)))))

+

+(defmethod expt :default [base pow] (Math/pow base pow))

+

+(defn abs "(abs n) is the absolute value of n" [n]

+  (cond

+   (not (number? n)) (throw (IllegalArgumentException.

+			     "abs requires a number"))

+   (neg? n) (- n)

+   :else n))

+

+(defmulti floor class)

+(defmethod floor ::integer [n] n)

+(defmethod floor java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_FLOOR) (toBigInteger)))

+(defmethod floor clojure.lang.Ratio [n]

+  (if (pos? n) (quot (. n numerator) (. n denominator))

+      (dec (quot (. n numerator) (. n denominator)))))

+(defmethod floor :default [n]

+  (Math/floor n))

+

+(defmulti ceil class)

+(defmethod ceil ::integer [n] n)

+(defmethod ceil java.math.BigDecimal [n] (.. n (setScale 0 BigDecimal/ROUND_CEILING) (toBigInteger)))

+(defmethod ceil clojure.lang.Ratio [n]

+  (if (pos? n) (inc (quot (. n numerator) (. n denominator)))

+      (quot (. n numerator) (. n denominator))))

+(defmethod ceil :default [n]

+  (Math/ceil n))

+

+(defmulti round class)

+(defmethod round ::integer [n] n)

+(defmethod round java.math.BigDecimal [n] (floor (+ n 0.5M)))

+(defmethod round clojure.lang.Ratio [n] (floor (+ n 1/2)))

+(defmethod round :default [n] (Math/round n))

+

+; mod is private because it is likely to be incorporated into Clojure's core

+(defn- mod [a b]

+  (cond

+   (or (not (integer? a)) (not (integer? b))) (throw (IllegalArgumentException.

+						      "mod requires two integers"))

+

+   (or (< a 0 b) (< b 0 a)) (+ (rem a b) b)

+   :else (rem a b)))

+

+(defn gcd "(gcd a b) returns the greatest common divisor of a and b" [a b]

+  (if (or (not (integer? a)) (not (integer? b)))

+    (throw (IllegalArgumentException. "gcd requires two integers"))  

+    (loop [a (abs a) b (abs b)]

+      (if (zero? b) a,

+	  (recur b (mod a b))))))

+

+; Length of integer in binary, used as helper function for sqrt.

+(defmulti #^{:private true} integer-length class)

+(defmethod integer-length java.lang.Integer [n]

+  (count (Integer/toBinaryString n)))

+(defmethod integer-length java.lang.Long [n]

+  (count (Integer/toBinaryString n)))

+(defmethod integer-length java.math.BigInteger [n]

+  (count (. n toString 2)))

+

+;; Produces the largest integer less than or equal to the square root of n

+;; Input n must be a non-negative integer

+(defn- integer-sqrt [n]

+  (cond

+   (> n 24)

+   (let [n-len (integer-length n)]

+     (loop [init-value (if (even? n-len)

+			 (bit-shift-left 1 (bit-shift-right n-len 1))

+			 (bit-shift-left 2 (bit-shift-right n-len 1)))]

+       (let [iterated-value (bit-shift-right (+ init-value (quot n init-value)) 1)]

+	 (if (>= iterated-value init-value)

+	   init-value

+	   (recur iterated-value)))))

+   (> n 15) 4

+   (> n  8) 3

+   (> n  3) 2

+   (> n  0) 1

+   (> n -1) 0))

+

+(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'.

+For example, (exact-integer-sqrt 15) is [3 6] because 15 = 3^2+6."

+  [n]

+  (if (or (not (integer? n)) (neg? n))

+    (throw (IllegalArgumentException. "exact-integer-sqrt requires a non-negative integer"))

+    (let [isqrt (integer-sqrt n),

+	  error (- n (* isqrt isqrt))]

+      [isqrt error])))

+

+(defmulti sqrt class)

+(defmethod sqrt ::integer [n]

+  (if (neg? n) Double/NaN

+      (let [isqrt (integer-sqrt n),

+	    error (- n (* isqrt isqrt))]

+	(if (zero? error) isqrt

+	    (Math/sqrt n)))))

+

+(defmethod sqrt clojure.lang.Ratio [n]

+  (if (neg? n) Double/NaN

+      (let [numerator (.numerator n),

+	    denominator (.denominator n),

+	    sqrtnum (sqrt numerator)]

+	(if (float? sqrtnum)

+	  (Math/sqrt n)

+	  (let [sqrtden (sqrt denominator)]

+	    (if (float? sqrtnum)

+	      (Math/sqrt n)

+	      (/ sqrtnum sqrtden)))))))

+

+(defmethod sqrt java.math.BigDecimal [n]

+  (if (neg? n) Double/NaN

+      (let [frac (rationalize n),

+	    sqrtfrac (sqrt frac)]

+	(if (ratio? sqrtfrac)

+	  (/ (BigDecimal. (.numerator sqrtfrac))

+	     (BigDecimal. (.denominator sqrtfrac)))

+	  sqrtfrac))))

+

+(defmethod sqrt :default [n]

+  (Math/sqrt n))
\ No newline at end of file
diff --git a/src/clojure/contrib/math/tests.clj b/src/clojure/contrib/math/tests.clj
new file mode 100644
index 00000000..19304924
--- /dev/null
+++ b/src/clojure/contrib/math/tests.clj
@@ -0,0 +1,96 @@
+(ns clojure.contrib.math.tests

+  (:use clojure.contrib.test-is

+	clojure.contrib.math))

+

+(deftest test-expt

+  (are (= _1 _2)

+      (expt 2 3) 8

+      (expt (expt 2 32) 2) (expt 2 64)

+      (expt 4/3 2) 16/9

+      (expt 2 -10) 1/1024

+      (expt 0.5M 2) 0.25M

+      (expt 5 4.2) (Math/pow 5 4.2)

+      (expt 5.3 4) (Math/pow 5.3 4)))

+

+(deftest test-abs

+  (are (= _1 _2)  

+      (abs -2) 2

+      (abs 0) 0

+      (abs 5) 5

+      (abs 123456789123456789) 123456789123456789

+      (abs -123456789123456789) 123456789123456789

+      (abs 5/3) 5/3

+      (abs -4/3) 4/3

+      (abs 4.3M) 4.3M

+      (abs -4.3M) 4.3M

+      (abs 2.8) 2.8

+      (abs -2.8) 2.8))

+

+(deftest test-gcd

+  (are (= _1 _2)

+      (gcd 4 3) 1

+      (gcd 24 12) 12

+      (gcd 24 27) 3))

+

+(deftest test-floor

+  (are (= _1 _2)

+      (floor 6) 6

+      (floor -6) -6

+      (floor 123456789123456789) 123456789123456789

+      (floor -123456789123456789) -123456789123456789

+      (floor 4/3) 1

+      (floor -4/3) -2

+      (floor 4.3M) 4

+      (floor -4.3M) -5

+      (floor 4.3) 4.0

+      (floor -4.3) -5.0))

+

+(deftest test-ceil

+  (are (= _1 _2)

+      (ceil 6) 6

+      (ceil -6) -6

+      (ceil 123456789123456789) 123456789123456789

+      (ceil -123456789123456789) -123456789123456789

+      (ceil 4/3) 2

+      (ceil -4/3) -1

+      (ceil 4.3M) 5

+      (ceil -4.3M) -4

+      (ceil 4.3) 5.0

+      (ceil -4.3) -4.0))

+

+(deftest test-round

+  (are (= _1 _2)

+      (round 6) 6

+      (round -6) -6

+      (round 123456789123456789) 123456789123456789

+      (round -123456789123456789) -123456789123456789

+      (round 4/3) 1

+      (round 5/3) 2

+      (round 5/2) 3

+      (round -4/3) -1

+      (round -5/3) -2

+      (round -5/2) -2

+      (round 4.3M) 4

+      (round 4.7M) 5

+      (round -4.3M) -4

+      (round -4.7M) -5

+      (round 4.5M) 5

+      (round -4.5M) -4

+      (round 4.3) 4

+      (round 4.7) 5

+      (round -4.3) -4

+      (round -4.7) -5

+      (round 4.5) 5

+      (round -4.5) -4))

+

+(deftest test-sqrt

+  (are (= _1 _2)

+      (sqrt 9) 3

+      (sqrt 16/9) 4/3

+      (sqrt 0.25M) 0.5M

+      (sqrt 2) (Math/sqrt 2)))

+

+(deftest test-exact-integer-sqrt

+  (are (= _1 _2)

+   (exact-integer-sqrt 15) [3 6]

+   (exact-integer-sqrt (inc (expt 2 64))) [(expt 2 32) 1]))