Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 249 lines (219 sloc) 8.739 kb
f433c7b @Engelberg Initial commit
Engelberg authored
1 ;;; math.clj: math functions that deal intelligently with the various
2 ;;; types in Clojure's numeric tower, as well as math functions
3 ;;; commonly found in Scheme implementations.
4
5 ;; by Mark Engelberg (mark.engelberg@gmail.com)
6 ;; May 21, 2011
7
8 (ns
9 ^{:author "Mark Engelberg",
10 :doc "Math functions that deal intelligently with the various
11 types in Clojure's numeric tower, as well as math functions
12 commonly found in Scheme implementations.
13
14 expt - (expt x y) is x to the yth power, returns an exact number
15 if the base is an exact number, and the power is an integer,
16 otherwise returns a double.
17 abs - (abs n) is the absolute value of n
18 gcd - (gcd m n) returns the greatest common divisor of m and n
19 lcm - (lcm m n) returns the least common multiple of m and n
20
21 When floor, ceil, and round are passed doubles, we just defer to
22 the corresponding functions in Java's Math library. Java's
23 behavior is somewhat strange (floor and ceil return doubles rather
24 than integers, and round on large doubles yields spurious results)
25 but it seems best to match Java's semantics. On exact numbers
26 (ratios and decimals), we can have cleaner semantics.
27
28 floor - (floor n) returns the greatest integer less than or equal to n.
29 If n is an exact number, floor returns an integer,
30 otherwise a double.
31 ceil - (ceil n) returns the least integer greater than or equal to n.
32 If n is an exact number, ceil returns an integer,
33 otherwise a double.
34 round - (round n) rounds to the nearest integer.
35 round always returns an integer. round rounds up for values
36 exactly in between two integers.
37
38
39 sqrt - Implements the sqrt behavior I'm accustomed to from PLT Scheme,
40 specifically, if the input is an exact number, and is a square
41 of an exact number, the output will be exact. The downside
42 is that for the common case (inexact square root), some extra
43 computation is done to look for an exact square root first.
44 So if you need blazingly fast square root performance, and you
45 know you're just going to need a double result, you're better
46 off calling java's Math/sqrt, or alternatively, you could just
47 convert your input to a double before calling this sqrt function.
48 If Clojure ever gets complex numbers, then this function will
49 need to be updated (so negative inputs yield complex outputs).
50 exact-integer-sqrt - Implements a math function from the R6RS Scheme
51 standard. (exact-integer-sqrt k) where k is a non-negative integer,
52 returns [s r] where k = s^2+r and k < (s+1)^2. In other words, it
53 returns the floor of the square root and the \"remainder\".
54 "}
55 clojure.math.numeric-tower)
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
56
57 ;; so this code works with both 1.2.x and 1.3.0:
58 (def ^{:private true} minus (first [-' -]))
59 (def ^{:private true} mult (first [*' *]))
60 (def ^{:private true} plus (first [+' +]))
61 (def ^{:private true} dec* (first [dec' dec]))
62 (def ^{:private true} inc* (first [inc' inc]))
f433c7b @Engelberg Initial commit
Engelberg authored
63
64 (defn- expt-int [base pow]
65 (loop [n pow, y (num 1), z base]
66 (let [t (even? n), n (quot n 2)]
67 (cond
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
68 t (recur n y (mult z z))
69 (zero? n) (mult z y)
70 :else (recur n (mult z y) (mult z z))))))
f433c7b @Engelberg Initial commit
Engelberg authored
71
72 (defn expt
73 "(expt base pow) is base to the pow power.
74 Returns an exact number if the base is an exact number and the power is an integer, otherwise returns a double."
75 [base pow]
76 (if (and (not (float? base)) (integer? pow))
77 (cond
78 (pos? pow) (expt-int base pow)
79 (zero? pow) 1
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
80 :else (/ 1 (expt-int base (minus pow))))
f433c7b @Engelberg Initial commit
Engelberg authored
81 (Math/pow base pow)))
82
83 (defn abs "(abs n) is the absolute value of n" [n]
84 (cond
85 (not (number? n)) (throw (IllegalArgumentException.
86 "abs requires a number"))
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
87 (neg? n) (minus n)
f433c7b @Engelberg Initial commit
Engelberg authored
88 :else n))
89
90 (defprotocol MathFunctions
91 (floor [n] "(floor n) returns the greatest integer less than or equal to n.
92 If n is an exact number, floor returns an integer, otherwise a double.")
93 (ceil [n] "(ceil n) returns the least integer greater than or equal to n.
94 If n is an exact number, ceil returns an integer, otherwise a double.")
95 (round [n] "(round n) rounds to the nearest integer.
96 round always returns an integer. Rounds up for values exactly in between two integers.")
97 (integer-length [n] "Length of integer in binary")
98 (sqrt [n] "Square root, but returns exact number if possible."))
99
100 (declare sqrt-integer)
101 (declare sqrt-ratio)
102 (declare sqrt-decimal)
103
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
104 ;; feature testing macro, based on suggestion from Chas Emerick:
105 (defmacro when-available
106 [sym & body]
107 (try
108 (when (resolve sym)
109 (list* 'do body))
110 (catch ClassNotFoundException _#)))
111
f433c7b @Engelberg Initial commit
Engelberg authored
112 (extend-type
113 Integer MathFunctions
114 (floor [n] n)
115 (ceil [n] n)
116 (round [n] n)
117 (integer-length [n] (- 32 (Integer/numberOfLeadingZeros n)))
118 (sqrt [n] (sqrt-integer n)))
119
120 (extend-type
121 Long MathFunctions
122 (floor [n] n)
123 (ceil [n] n)
124 (round [n] n)
125 (integer-length [n] (- 64 (Long/numberOfLeadingZeros n)))
126 (sqrt [n] (sqrt-integer n)))
127
128 (extend-type
129 java.math.BigInteger MathFunctions
130 (floor [n] n)
131 (ceil [n] n)
132 (round [n] n)
133 (integer-length [n] (.bitLength n))
134 (sqrt [n] (sqrt-integer n)))
135
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
136 (when-available
137 clojure.lang.BigInt
138 (extend-type
139 clojure.lang.BigInt MathFunctions
140 (floor [n] n)
141 (ceil [n] n)
142 (round [n] n)
143 (integer-length [n] (.bitLength n))
144 (sqrt [n] (sqrt-integer n))))
f433c7b @Engelberg Initial commit
Engelberg authored
145
146 (extend-type
147 java.math.BigDecimal MathFunctions
148 (floor [n] (bigint (.setScale n 0 BigDecimal/ROUND_FLOOR)))
149 (ceil [n] (bigint (.setScale n 0 BigDecimal/ROUND_CEILING)))
150 (round [n] (floor (+ n 0.5M)))
151 (sqrt [n] (sqrt-decimal n)))
152
153 (extend-type
154 clojure.lang.Ratio MathFunctions
155 (floor [n]
156 (if (pos? n) (quot (. n numerator) (. n denominator))
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
157 (dec* (quot (. n numerator) (. n denominator)))))
f433c7b @Engelberg Initial commit
Engelberg authored
158 (ceil [n]
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
159 (if (pos? n) (inc* (quot (. n numerator) (. n denominator)))
f433c7b @Engelberg Initial commit
Engelberg authored
160 (quot (. n numerator) (. n denominator))))
161 (round [n] (floor (+ n 1/2)))
162 (sqrt [n] (sqrt-ratio n)))
163
164 (extend-type
165 Double MathFunctions
166 (floor [n] (Math/floor n))
167 (ceil [n] (Math/ceil n))
168 (round [n] (Math/round n)) ;(round (bigdec n)))
169 (sqrt [n] (Math/sqrt n)))
170
171 (extend-type
172 Float MathFunctions
173 (floor [n] (Math/floor n))
174 (ceil [n] (Math/ceil n))
175 (round [n] (Math/round n)) ;(round (bigdec n)))
176 (sqrt [n] (Math/sqrt n)))
177
178 (defn gcd "(gcd a b) returns the greatest common divisor of a and b" [a b]
179 (if (or (not (integer? a)) (not (integer? b)))
180 (throw (IllegalArgumentException. "gcd requires two integers"))
181 (loop [a (abs a) b (abs b)]
182 (if (zero? b) a,
183 (recur b (mod a b))))))
184
185 (defn lcm
186 "(lcm a b) returns the least common multiple of a and b"
187 [a b]
188 (when (or (not (integer? a)) (not (integer? b)))
189 (throw (IllegalArgumentException. "lcm requires two integers")))
190 (cond (zero? a) 0
191 (zero? b) 0
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
192 :else (abs (mult b (quot a (gcd a b))))))
f433c7b @Engelberg Initial commit
Engelberg authored
193
194 ;; Produces the largest integer less than or equal to the square root of n
195 ;; Input n must be a non-negative integer
196 (defn- integer-sqrt [n]
197 (cond
198 (> n 24)
199 (let [n-len (integer-length n)]
200 (loop [init-value (if (even? n-len)
201 (expt 2 (quot n-len 2))
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
202 (expt 2 (inc* (quot n-len 2))))]
203 (let [iterated-value (quot (plus init-value (quot n init-value)) 2)]
f433c7b @Engelberg Initial commit
Engelberg authored
204 (if (>= iterated-value init-value)
205 init-value
206 (recur iterated-value)))))
207 (> n 15) 4
208 (> n 8) 3
209 (> n 3) 2
210 (> n 0) 1
211 (> n -1) 0))
212
213 (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'.
214 For example, (exact-integer-sqrt 15) is [3 6] because 15 = 3^2+6."
215 [n]
216 (if (or (not (integer? n)) (neg? n))
217 (throw (IllegalArgumentException. "exact-integer-sqrt requires a non-negative integer"))
218 (let [isqrt (integer-sqrt n),
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
219 error (minus n (mult isqrt isqrt))]
f433c7b @Engelberg Initial commit
Engelberg authored
220 [isqrt error])))
221
222 (defn- sqrt-integer [n]
223 (if (neg? n) Double/NaN
224 (let [isqrt (integer-sqrt n),
6a961d5 @seancorfield Clojure 1.2.x compatibility
seancorfield authored
225 error (minus n (mult isqrt isqrt))]
f433c7b @Engelberg Initial commit
Engelberg authored
226 (if (zero? error) isqrt
227 (Math/sqrt n)))))
228
229 (defn- sqrt-ratio [n]
230 (if (neg? n) Double/NaN
231 (let [numerator (.numerator n),
232 denominator (.denominator n),
233 sqrtnum (sqrt numerator)]
234 (if (float? sqrtnum)
235 (Math/sqrt n)
236 (let [sqrtden (sqrt denominator)]
237 (if (float? sqrtnum)
238 (Math/sqrt n)
239 (/ sqrtnum sqrtden)))))))
240
241 (defn- sqrt-decimal [n]
242 (if (neg? n) Double/NaN
243 (let [frac (rationalize n),
244 sqrtfrac (sqrt frac)]
245 (if (ratio? sqrtfrac)
246 (/ (BigDecimal. (.numerator sqrtfrac))
247 (BigDecimal. (.denominator sqrtfrac)))
248 sqrtfrac))))
Something went wrong with that request. Please try again.