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