Last active
October 5, 2018 23:11
-
-
Save adamdavislee/d7681282fafefc17929ece0d291c985f to your computer and use it in GitHub Desktop.
Translation of the Number Theory Namespace of Racket's Math Library: https://docs.racket-lang.org/math/number-theory.html
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
(ns adamdavislee.math.number-theory | |
"This port closely mimics the naming conventions and coding style of the racket source code. | |
However it does not: | |
Define inline versions of mod+, mod*, etcetera. | |
TODO: | |
Add pre/post conditions and assertions with type checking etc. | |
Improve docstrings. | |
Finish marking internal functions as private. | |
Use internal memoization for some of the functions which racket uses manual bit manipulation for." | |
(:require [clojure.test | |
:refer [is | |
are | |
with-test | |
successful? | |
run-tests]])) | |
(do | |
(defmacro | |
^:private | |
defval | |
"like def, but returns value instead of var" | |
[sym value] | |
`(do (def ~sym ~value) ~sym)) | |
(with-test | |
#'adamdavislee.math.number-theory/defval | |
(is | |
(and (keyword? (defval a :a)) | |
(var? (def a :a)) | |
(do (defval b :b) | |
(= b :b)))))) | |
(with-test | |
(defn | |
assert | |
[value predicate & [error-message]] | |
(if (predicate value) value error-message)) | |
(is (nil? (assert "" int?))) | |
(is (= "oops!" (assert "" int? "oops!"))) | |
(is (= "" (assert "" string?)))) | |
(with-test | |
(defn intvalue? [n] (== n (int n))) | |
(are | |
[a b] | |
(= (intvalue? a) b) | |
1.2 | |
false | |
1.0 | |
true | |
-3.0 | |
true | |
-3.3 | |
false | |
-2 | |
true)) | |
(with-test | |
(defn- abs "absolute value" [a] (if (neg? a) (- a) a)) | |
(is (= (abs -1) 1)) | |
(is (= (abs 1) 1))) | |
(with-test | |
(defn- | |
pow | |
([] 1) | |
([n] n) | |
([a b] | |
(if | |
(intvalue? b) | |
((if (neg? b) / *) (apply * (repeat (abs b) a))) | |
(Math/pow a b))) | |
([a b & numbers] | |
(let [numbers (conj numbers b a)] (reduce pow numbers)))) | |
(are | |
[a b] | |
(= (double (apply pow a)) (double b)) | |
[] | |
1 | |
[2] | |
2 | |
[2 2] | |
4 | |
[2 2 2] | |
16 | |
[2 3 4] | |
(pow (pow 2 3) 4) | |
[3 -1] | |
1/3 | |
[4 -2] | |
1/16)) | |
(with-test (defn- rt ([a b] (pow a (/ b)))) (is (= (rt 8 3) 2.0))) | |
(with-test | |
(def | |
one? | |
"imo it's reasonable to introduce vocabulary for 0 1 and 2" | |
(partial = 1)) | |
(are [a b] (= (apply one? a) b) [1] true [2] false)) | |
(with-test (def duple (partial * 2)) (is (= 6 (duple 3)))) | |
(with-test (def two? (partial = 2)) (is (not (two? 1))) (is (two? 2))) | |
(with-test (defn- half [n] (/ n 2)) (is (= (half 4) 2))) | |
(with-test | |
(defn- sqrt [n] (Math/sqrt n)) | |
(is (= (sqrt 1) (double 1)))) | |
(with-test (defn- sq [n] (pow n 2)) (is (= (sq 3) 9))) | |
(with-test | |
(defn- | |
gcd | |
"uses the euclidean algorithm, see https://brilliant.org/wiki/euclidean-algorithm" | |
([] 0) | |
([a] a) | |
([a b] (if (zero? b) (abs a) (recur b (rem a b)))) | |
([a b & numbers] (gcd a (apply gcd (conj numbers b))))) | |
(are | |
[a b] | |
(= (apply gcd a) b) | |
[] | |
0 | |
[3] | |
3 | |
[2 4] | |
2 | |
[3 4] | |
1 | |
[-2 4] | |
2 | |
[3 4 5] | |
1)) | |
(with-test | |
(defn- | |
inrange | |
"range with an inclusive end and default beginning of 1" | |
([] (drop 1 (range))) | |
([end] (range 1 (inc end))) | |
([start end] (range start (inc end))) | |
([start end step] (range start (inc end) step))) | |
(are | |
[a b] | |
(= (apply inrange a) b) | |
[8] | |
[1 2 3 4 5 6 7 8] | |
[2 4] | |
[2 3 4])) | |
(with-test | |
(defn- andmap [& args] (every? identity (apply map args))) | |
(are | |
[a b] | |
(= (apply andmap a) b) | |
[pos? []] | |
true | |
[neg? [-2 -1]] | |
true | |
[neg? [-2 -1 0 1 2]] | |
false)) | |
(def ^:private fact-table-size 171) | |
(def | |
fact-table | |
(vec | |
(reverse | |
(reduce | |
(fn [numbers n] (cons (* n (first numbers)) numbers)) | |
[1N] | |
(inrange fact-table-size))))) | |
(def ^:private simple-cutoff 244) | |
(with-test | |
(defn- | |
factorial-simple | |
[n] | |
(if | |
(< n fact-table-size) | |
(fact-table n) | |
(* n (factorial-simple (dec n))))) | |
(is | |
(andmap | |
(fn [n] (= (factorial-simple n) (apply * (inrange n)))) | |
(inrange 8)))) | |
(with-test | |
(defn divides? [a b] (if (zero? a) false (= 0 (rem b a)))) | |
(are [a b] (divides? a b) 2 4 2 12 2 0) | |
(are [a b] (not (divides? a b)) 0 1 2 3 2 13)) | |
(with-test | |
(defn coprime? [a & numbers] (one? (apply gcd (conj numbers a)))) | |
(are | |
[a b] | |
(= (apply coprime? a) b) | |
[1] | |
true | |
[2] | |
false | |
[2 9] | |
true | |
[2 4 8] | |
false | |
[2 4 9] | |
true) | |
(is (coprime? (* 3 7) (* 5 19))) | |
(is (not (coprime? (* 3 7 5) (* 5 19 2))))) | |
(with-test | |
(defn | |
pairwise-coprime? | |
[a & numbers] | |
(or | |
(nil? numbers) | |
(and | |
(if | |
(every? identity (map (fn [b] (coprime? a b)) numbers)) | |
true | |
false) | |
(if (apply pairwise-coprime? numbers) true false)))) | |
(are | |
[a b] | |
(= (apply pairwise-coprime? a) b) | |
[1] | |
true | |
[1 2] | |
true | |
[1 2 4] | |
false | |
[1 2 -4] | |
false) | |
(are [a] (not (apply pairwise-coprime? a)) [10 7 33 14] [6 10 15]) | |
(is (pairwise-coprime? 10 7 33 13))) | |
(with-test | |
(defn | |
bezout | |
"uses the extended euclidean algorithm, see: https://brilliant.org/wiki/extended-euclidean-algorithm" | |
([a] 1) | |
([a b] | |
(let | |
[loop | |
(fn | |
[a b ua va ub vb] | |
(let | |
[q (quot a b) r (rem a b)] | |
(if | |
(= r 0) | |
(list ub vb) | |
(recur b r ub vb (- ua (* q ub)) (- va (* q vb)))))) | |
start | |
(fn [a b] (if (> a b) (loop a b 1 0 0 1) (loop b a 0 1 1 0)))] | |
(cond | |
(and (pos? a) (pos? b)) | |
(start a b) | |
(and (neg? a) (neg? b)) | |
(mapv - (start (- a) (- b))) | |
(and (neg? a) (pos? b)) | |
(let | |
[k (+ (quot (- a) b) 1) [u v] (start (+ a (* k b)) b)] | |
[u (+ (* u k) v)]) | |
(and (pos? a) (neg? b)) | |
(let | |
[k (+ (quot (- b) a) 1) [u v] (start a (+ (* k a) b))] | |
(list (+ u (* k v)) v))))) | |
([a b & numbers] | |
(let | |
[numbers | |
(conj numbers b a) | |
[s t] | |
(bezout (apply gcd (rest numbers)) (first numbers))] | |
(conj (map (partial * s) (apply bezout (rest numbers))) t)))) | |
(are | |
[a b] | |
(and (= (apply bezout a) b) (= (apply gcd a) (apply + (map * a b)))) | |
[4 2] | |
[0 1] | |
[2 4] | |
[1 0] | |
[2 -4] | |
[1 0] | |
[-2 4] | |
[1 1] | |
[-2 -4] | |
[-1 0] | |
[3 4] | |
[-1 1] | |
[2 4 8] | |
[1 0 0] | |
[2 4 7] | |
[0 2 -1]) | |
(are | |
[a] | |
(= (apply gcd a) (apply + (map * (apply bezout a) a))) | |
[12 20] | |
[12 20] | |
[20 16] | |
[12 20 16])) | |
(with-test | |
(defn modular-inverse [a n] (mod (first (bezout a n)) n)) | |
(are [a b] (= (apply modular-inverse a) b) [3 7] 5) | |
(andmap | |
(fn | |
[n] | |
(let | |
[m (and (coprime? n 20) (modular-inverse n 20))] | |
(if m (= (rem (* n m) 20) 1) true))) | |
(inrange 20))) | |
(with-test | |
(defn | |
solve-chinese | |
[as ns] | |
(let | |
[n | |
(apply * ns) | |
cs | |
(map (fn [ni] (quot n ni)) ns) | |
ds | |
(map modular-inverse cs ns)] | |
(mod (apply + (map * as cs ds)) n))) | |
(are [a b] (= (apply solve-chinese a) b) [[2 3 2] [3 5 7]] 23)) | |
(with-test | |
(defn | |
mediant | |
[a b] | |
(/ | |
(+ (numerator a) (numerator b)) | |
(+ (denominator a) (denominator b)))) | |
(are [a b] (= (apply mediant a) b) [1/2 5/6] 3/4 [-1/2 -2/3] -3/5)) | |
(with-test | |
(defn | |
perfect-square? | |
[a] | |
(let [a-sqrt (int (sqrt a))] (if (= (pow a-sqrt 2) a) a-sqrt false))) | |
(are [a b] (= (apply perfect-square? a) b) [4] 2 [5] false)) | |
(with-test | |
(defn triangle-number [n] (quot (* n (inc n)) 2)) | |
(are [a b] (= (apply triangle-number a) b) [1] 1 [2] 3) | |
(is (= (mapv triangle-number (inrange 0 5)) [0 1 3 6 10 15]))) | |
(with-test | |
(defn square-number? [n] (and (perfect-square? n) true)) | |
(is (andmap square-number? [0 1 4 9 16 25])) | |
(is (andmap (complement square-number?) [2 5 10 17 26]))) | |
(with-test | |
(defn pentagonal-number [n] (quot (* n (- (* 3 n) 1)) 2)) | |
(is (= (mapv pentagonal-number [0 1 2 3 4 5]) [0 1 5 12 22 35]))) | |
(with-test | |
(defn hexagonal-number [n] (* n (- (* 2 n) 1))) | |
(is (= (mapv hexagonal-number [0 1 2 3 4 5]) [0 1 6 15 28 45]))) | |
(with-test | |
(defn heptagonal-number [n] (quot (* n (- (* 5 n) 3)) 2)) | |
(is (= (mapv heptagonal-number [0 1 2 3 4 5]) [0 1 7 18 34 55]))) | |
(with-test | |
(defn octagonal-number [n] (* n (- (* 3 n) 2))) | |
(is (= (mapv octagonal-number [0 1 2 3 4 5]) [0 1 8 21 40 65]))) | |
(with-test | |
(defn | |
quadratic-solutions | |
"return list of solutions to a x^2 + b x + c = 0" | |
[a b c] | |
(let | |
[d (- (* b b) (* 4 a c))] | |
(cond | |
(< d 0) | |
[] | |
(= d 0) | |
[(/ b (* -2 a))] | |
:else | |
(let | |
[sqrt-d (sqrt d)] | |
[(/ (- (- b) sqrt-d) (* 2 a)) (/ (+ (- b) sqrt-d) (* 2 a))])))) | |
(are | |
[a] | |
a | |
(= (quadratic-solutions 1 0 -4) [-2.0 2.0]) | |
(= (quadratic-solutions 1 0 4) []) | |
(= (quadratic-solutions 1 0 0) [0]))) | |
(with-test | |
(defn exact-integer? [n] (= (double n) (double (int n)))) | |
(is (exact-integer? 1.0)) | |
(is (not (exact-integer? 1.1)))) | |
(with-test | |
(defn | |
quadratic-integer-solutions | |
[a b c] | |
"return list of integer solutions to a x^2 + b x + c = 0" | |
(filter exact-integer? (quadratic-solutions a b c))) | |
(are | |
[a b] | |
(= (apply quadratic-integer-solutions a) b) | |
[1 0 -4] | |
[-2.0 2.0] | |
[1 0 4] | |
[] | |
[1 0 0] | |
[0])) | |
(with-test | |
(defn | |
natural? | |
[n] | |
(every? | |
identity | |
((juxt | |
exact-integer? | |
(fn [n] (some identity ((juxt pos? zero?) n)))) | |
n))) | |
(are [a] (natural? a) 0 1 2 3 4) | |
(are [a] (not (natural? a)) -1 -2 -3 -4)) | |
(with-test | |
(defn | |
quadratic-natural-solutions | |
"return list of nonnegative-integer solutions to a x^2 + b x + c = 0" | |
[a b c] | |
(filter natural? (quadratic-solutions a b c))) | |
(are | |
[a b] | |
(= (apply quadratic-natural-solutions a) b) | |
[1 0 -4] | |
[2.0] | |
[1 0 4] | |
[] | |
[1 0 0] | |
[0])) | |
(with-test | |
(defn | |
tangent-number | |
"the nth tangent number, http://mathworld.wolfram.com/TangentNumber.html" | |
[n] | |
(def T (atom (assoc (vec (repeat (+ n 2) 0)) 1 1))) | |
(doseq | |
[k (range (inc n))] | |
(do | |
(doseq | |
[i (range (inc k))] | |
(swap! | |
T | |
(fn* | |
[p1__21131#] | |
(assoc p1__21131# i (* (inc i) (p1__21131# (inc i))))))) | |
(swap! T (fn* [p1__21132#] (assoc p1__21132# k 0))) | |
(doseq | |
[i (range (inc k) 1 -1)] | |
(swap! | |
T | |
(fn* | |
[p1__21133#] | |
(assoc p1__21133# i (+ (p1__21133# i) (p1__21133# (- i 2))))))))) | |
(first (clojure.core/deref T))) | |
(is | |
(= | |
(mapv tangent-number [1 3 5 7 9 11 13]) | |
[1 2 16 272 7936 353792 22368256])) | |
(is (= (mapv tangent-number [0 2 4 6 8 10]) [0 0 0 0 0 0]))) | |
(with-test | |
(defn | |
farey-sequence | |
[n] | |
(loop | |
[a 1 b 1 c (dec n) d n fs []] | |
(let | |
[fs (cons (/ a b) fs)] | |
(if | |
(pos? a) | |
(let | |
[k (quot (+ n b) d)] | |
(recur c d (- (* k c) a) (- (* k d) b) fs)) | |
fs)))) | |
(is | |
(= | |
(farey-sequence 5) | |
(quote (0 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 1))))) | |
(with-test | |
(defn | |
eulerian-number | |
"uses standard recurrence: <n,k> = (k+1) <n-1,k> + (n-k) <n-1,k-1>" | |
[n k] | |
(if | |
(zero? k) | |
1 | |
(do | |
(def E (atom (vec (repeat (max (inc n) (inc k)) 0)))) | |
(swap! E (fn* [p1__21134#] (assoc p1__21134# 0 1))) | |
(doseq | |
[i (inrange n) j (range (dec i) 0 -1)] | |
(swap! | |
E | |
(fn* | |
[p1__21135#] | |
(assoc | |
p1__21135# | |
j | |
(+ | |
(* (inc j) ((clojure.core/deref E) j)) | |
(* (- i j) ((clojure.core/deref E) (dec j)))))))) | |
((clojure.core/deref E) k)))) | |
(is | |
(= (mapv (partial eulerian-number 5) [0 1 2 3 4]) [1 26 66 26 1]))) | |
(with-test | |
(defn | |
binomial | |
[n k] | |
(if | |
(or (zero? k) (every? one? [n k])) | |
1 | |
(loop | |
[n n k k] | |
(cond | |
(zero? k) | |
1 | |
(one? k) | |
n | |
(> k n) | |
0 | |
(two? k) | |
(half (* n (dec n))) | |
(> k (half n)) | |
(recur n (- n k)) | |
:else | |
(* | |
(+ n (- k) 1) | |
(reduce | |
(fn [product i] (* product (/ (+ n (- k) i) i))) | |
1 | |
(inrange 2 k))))))) | |
(are | |
[a b] | |
(= (apply binomial a) b) | |
[10 3] | |
120 | |
[10 11] | |
0 | |
[10 0] | |
1 | |
[10 10] | |
1 | |
[10 1] | |
10 | |
[10 9] | |
10)) | |
(with-test | |
(defn | |
triangle-number? | |
[n] | |
(not (nil? (quadratic-natural-solutions 1/2 1/2 (- n))))) | |
(is (= (andmap triangle-number? (quote (0 1 3 6 10 15)))))) | |
(with-test | |
(defn | |
pentagonal-number? | |
[n] | |
(not (nil? (quadratic-natural-solutions 3/2 -1/2 (- n))))) | |
(is (= (andmap pentagonal-number? (quote (0 1 5 12 22 35)))))) | |
(with-test | |
(defn | |
hexagonal-number? | |
[n] | |
(not (nil? (quadratic-natural-solutions 2 -1 (- n))))) | |
(is (= (andmap hexagonal-number? (quote (0 1 6 15 28 45)))))) | |
(with-test | |
(defn | |
heptagonal-number? | |
[n] | |
(not (nil? (quadratic-natural-solutions 5/2 -3/2 (- n))))) | |
(is (= (andmap heptagonal-number? (quote (0 1 7 18 34 55)))))) | |
(with-test | |
(defn | |
octagonal-number? | |
[n] | |
(not (nil? (quadratic-natural-solutions 3 -2 (- n))))) | |
(is (= (andmap octagonal-number? (quote (0 1 8 21 40 65)))))) | |
(with-test | |
(defn | |
make-fibonacci | |
[c d] | |
(fn | |
[n] | |
(loop | |
[a d b c p 0 q 1 count n] | |
(cond | |
(zero? count) | |
b | |
(even? count) | |
(recur | |
a | |
b | |
(+ (* p p) (* q q)) | |
(+ (* 2 p q) (* q q)) | |
(quot count 2)) | |
:else | |
(recur | |
(+ (* b q) (* a q) (* a p)) | |
(+ (* b p) (* a q)) | |
p | |
q | |
(dec count)))))) | |
(is (= (mapv (make-fibonacci 2 1) (range 8)) [2 1 3 4 7 11 18 29]))) | |
(with-test | |
(def fibonacci (make-fibonacci 0 1)) | |
(= (mapv fibonacci (range 8)) [0 1 1 2 3 5 8 13])) | |
(with-test | |
(defn | |
make-modular-fibonacci | |
[c d] | |
(fn | |
[n modulus] | |
(loop | |
[a d b c p 0 q 1 count n] | |
(cond | |
(zero? count) | |
(mod b modulus) | |
(even? count) | |
(recur | |
a | |
b | |
(mod (+ (* p p) (* q q)) modulus) | |
(mod (+ (* 2 p q) (* q q)) modulus) | |
(quot count 2)) | |
:else | |
(recur | |
(mod (+ (* b q) (* a q) (* a p)) modulus) | |
(mod (+ (* b p) (* a q)) modulus) | |
p | |
q | |
(dec count)))))) | |
(is | |
(every? | |
identity | |
(for | |
[a (range -5 6) b (range -5 6) modulus (range 1 8)] | |
(= | |
(mapv | |
(fn [n] ((make-modular-fibonacci a b) n modulus)) | |
(range 20)) | |
(mapv (fn [n] (mod ((make-fibonacci a b) n) modulus)) (range 20))))))) | |
(with-test | |
(def modular-fibonacci (make-modular-fibonacci 0 1)) | |
(is | |
(= | |
(mapv (fn [n] (modular-fibonacci n 3)) (range 8)) | |
[0 1 1 2 0 2 2 1]))) | |
(with-test | |
(defn | |
factorial | |
[n] | |
(cond | |
(or (zero? n) (one? n)) | |
1 | |
(< n simple-cutoff) | |
(factorial-simple n) | |
:else | |
((fn | |
loop | |
[n a] | |
(if | |
(<= (- n a) 0) | |
n | |
(* (loop n (duple a)) (loop (- n a) (duple a))))) | |
n | |
1N))) | |
(defn | |
test-factorial | |
[n] | |
(if (zero? n) 1 (* (bigint n) (test-factorial (dec n))))) | |
(are | |
[a b] | |
(= a b) | |
(map factorial [0 1 2 3 4 5]) | |
[1 1 2 6 24 120] | |
(factorial (inc fact-table-size)) | |
(test-factorial (inc fact-table-size)) | |
(factorial (inc simple-cutoff)) | |
(test-factorial (inc simple-cutoff)))) | |
(with-test | |
(defn | |
permutations | |
[n k] | |
(cond | |
(zero? k) | |
1 | |
(< n k) | |
0 | |
:else | |
(/ (factorial n) (factorial (- n k))))) | |
10 | |
3 | |
(are | |
[a b] | |
(= (apply permutations a) b) | |
[10 3] | |
720 | |
[10 0] | |
1 | |
[10 10] | |
3628800 | |
[0 0] | |
1)) | |
(with-test | |
(defn | |
multinomial | |
[n ks] | |
(if | |
(not= n (apply + ks)) | |
0 | |
(apply / (factorial n) (map factorial ks)))) | |
(are | |
[a b] | |
(= (apply multinomial a) b) | |
[20 [3 4 5 8]] | |
3491888400 | |
[0 []] | |
1 | |
[4 [1 1]] | |
0)) | |
(with-test | |
(defn | |
modular-expt | |
[a b n] | |
((fn | |
loop | |
[a b] | |
(cond | |
(<= b 1) | |
(if (zero? b) (mod 1 n) (mod a n)) | |
(even? b) | |
(mod (sq (loop a (quot b 2))) n) | |
:else | |
(mod (* a (loop a (dec b))) n))) | |
a | |
b)) | |
(is (= (mod (pow -6N 523N) 19) (modular-expt -6N 523N 19)))) | |
(with-test | |
(defn- | |
modular-const | |
[n a] | |
(if | |
(integer? a) | |
(mod a n) | |
(mod (* (numerator a) (modular-inverse n (denominator a))) n))) | |
(are | |
[a b] | |
(= (apply modular-const a) b) | |
[1 2] | |
0 | |
[3 1] | |
1 | |
[2 1] | |
1 | |
[1 2] | |
0 | |
[3 2/3] | |
2 | |
[3 1/3] | |
1)) | |
(def ^:private ^:dynamic current-modulus 1) | |
(do | |
(defmacro | |
with-modulus | |
"binds the current-modulus atom to n while executing body" | |
[n & body] | |
`(binding | |
[current-modulus ~n] | |
(do ~@body))) | |
(with-test | |
#'adamdavislee.math.number-theory/with-modulus | |
(let [prev-modulus current-modulus | |
n (rand-int 1024)] | |
(is | |
(= n | |
(with-modulus | |
n | |
current-modulus))) | |
(is | |
(= prev-modulus | |
current-modulus))))) | |
(with-test | |
(defn | |
mod+ | |
([] 0) | |
([a] (mod a current-modulus)) | |
([a b] (mod (+ a b) current-modulus)) | |
([a b & cs] | |
(reduce | |
(fn [c d] (mod (+ c d) current-modulus)) | |
(mod (+ a b) current-modulus) | |
cs))) | |
(with-modulus | |
3 | |
(are [a b] (= (mod+ a) b) 92 2 79 1 109 1 106 1) | |
(are | |
[a b] | |
(= (apply mod+ a) b) | |
[117 11 4] | |
0 | |
[109 74 94] | |
1 | |
[49 66 65] | |
0 | |
[62 112 77] | |
2))) | |
(with-test | |
(defn | |
mod- | |
([a] (mod (- a) current-modulus)) | |
([a b] (mod (- a b) current-modulus)) | |
([a b & cs] | |
(reduce | |
(fn [c d] (mod (- c d) current-modulus)) | |
(mod (- a b) current-modulus) | |
cs))) | |
(with-modulus | |
3 | |
(are | |
[a b] | |
(= (apply mod- a) b) | |
[46 109 15] | |
0 | |
[62 102 113] | |
0 | |
[58 92 49] | |
1 | |
[118 54 98] | |
2))) | |
(with-test | |
(defn | |
mod* | |
([] (mod 1 current-modulus)) | |
([a] (mod a current-modulus)) | |
([a b] (mod (* a b) current-modulus)) | |
([a b & cs] | |
(reduce | |
(fn [c d] (mod (* c d) current-modulus)) | |
(mod (* a b) current-modulus) | |
cs))) | |
(with-modulus | |
3 | |
(let | |
[ins | |
(vec | |
(repeatedly | |
4 | |
(fn [] (vec (repeatedly 3 (fn [] (rand-int 128))))))) | |
outs | |
(mapv (fn [a] (apply mod* a)) ins)] | |
[ins outs]) | |
(are | |
[a b] | |
(= (apply mod* a) b) | |
[122 7 68] | |
1 | |
[75 63 98] | |
0 | |
[77 94 16] | |
2 | |
[49 18 9] | |
0))) | |
(with-test | |
(defn | |
mod-div | |
([a] (modular-inverse a current-modulus)) | |
([a b] | |
(mod (* a (modular-inverse b current-modulus)) current-modulus)) | |
([a b & cs] (mod-div a (apply mod* b cs)))) | |
(with-modulus | |
3 | |
(are | |
[a b] | |
(= (apply mod-div a) b) | |
[24971 28175 57626] | |
2 | |
[47928 26377 64096] | |
0 | |
[13786 27692 6754] | |
2 | |
[63868 50846 51536] | |
1))) | |
(with-test | |
(defn modsqr [a] (mod (sq a) current-modulus)) | |
(with-modulus | |
3 | |
(are | |
[a b] | |
(= (modsqr a) b) | |
(modsqr 3) | |
0 | |
(modsqr 4) | |
1 | |
(modsqr 5) | |
1 | |
(modsqr 6) | |
0))) | |
(with-test | |
(defn | |
modexpt | |
[a b] | |
(if | |
(< b 0) | |
(modular-expt | |
(modular-inverse a current-modulus) | |
(- b) | |
current-modulus) | |
(modular-expt a b current-modulus))) | |
(with-modulus 3 (are [a b] (= (apply modexpt a) b) [2 3] 2 [4 7] 1))) | |
(with-test | |
(defn | |
mod2 | |
[a] | |
(if | |
(integer? a) | |
(mod a current-modulus) | |
(mod-div (numerator a) (denominator a)))) | |
(is (zero? (with-modulus 7 (mod2 (* 218 7))))) | |
(are [a b] (= (with-modulus 7 (mod2 a)) b) (* 218 7) 0 3/2 5)) | |
(with-test | |
(defn- | |
def-mod-op | |
[op] | |
(fn [& numbers] (apply op (map mod2 numbers)))) | |
(is (with-modulus 3 ((def-mod-op =) 0 3 6 9)))) | |
(with-test | |
(def mod= (def-mod-op =)) | |
(with-modulus 3 (is (mod= 3 6 9 12)) (is (not (mod= 3 6 9 13))))) | |
(with-test | |
(def mod< (def-mod-op <)) | |
(with-modulus 3 (is (mod< 3 4 5)) (is (not (mod< 3 4 5 6))))) | |
(with-test | |
(def mod> (def-mod-op >)) | |
(with-modulus 3 (is (mod> 5 4 3)) (is (not (mod> 6 5 4 3))))) | |
(with-test | |
(def mod<= (def-mod-op <=)) | |
(with-modulus 3 (is (mod<= 3 3 4 4)) (is (not (mod<= 4 4 3 3))))) | |
(with-test | |
(def mod>= (def-mod-op >=)) | |
(with-modulus 3 (is (mod>= 4 4 3 3)) (is (not (mod>= 3 3 4 4))))) | |
(with-test | |
(defn | |
bernoulli | |
"based on this recursive definition: https://en.wikipedia.org/wiki/Bernoulli_number#Recursive_definition" | |
[n] | |
(- | |
(if (zero? n) 1 0) | |
(apply | |
+ | |
(mapv | |
(fn [k] (* (binomial n k) (/ (bernoulli k) (+ n (- k) 1)))) | |
(inrange 0 (dec n)))))) | |
(is | |
(= | |
(mapv bernoulli (inrange 0 16)) | |
[1 | |
-1/2 | |
1/6 | |
0 | |
-1/30 | |
0 | |
1/42 | |
0 | |
-1/30 | |
0 | |
5/66 | |
0 | |
-691/2730 | |
0 | |
7/6 | |
0 | |
-3617/510]))) | |
(with-test | |
(defn | |
partitions | |
[n] | |
(if | |
(zero? n) | |
1 | |
(+ | |
(loop | |
[k 1 m (dec n) s 0] | |
(if | |
(neg? m) | |
s | |
(recur | |
(inc k) | |
(- m (+ (* 3 k) 1)) | |
(if (odd? k) (+ s (partitions m)) (- s (partitions m)))))) | |
(loop | |
[k -1 m (- n 2) s 0] | |
(if | |
(neg? m) | |
s | |
(recur | |
(dec k) | |
(+ m (* 3 k) -2) | |
(if (odd? k) (+ s (partitions m)) (- s (partitions m))))))))) | |
(is | |
(= | |
(map partitions [0 1 2 3 4 5 6 7 8 9 10]) | |
[1 1 2 3 5 7 11 15 22 30 42]))) | |
(with-test | |
(defn- | |
nth-prime | |
[n] | |
(loop | |
[numbers (drop 2 (inrange)) primes [2] prime 2] | |
(if | |
(= n (count primes)) | |
(last primes) | |
(let | |
[numbers | |
(filter | |
(fn [a] (or (= prime a) (not (divides? prime a)))) | |
numbers) | |
primes | |
(conj primes (first numbers)) | |
prime | |
(first numbers) | |
numbers | |
(rest numbers)] | |
(recur numbers primes prime))))) | |
(is (= (map nth-prime (inrange 1 8)) [2 3 5 7 11 13 17 19]))) | |
(with-test | |
(defn | |
prime? | |
[n] | |
(cond | |
(one? n) | |
false | |
(<= n 3) | |
true | |
(or (divides? 2 n) (divides? 3 n)) | |
false | |
:else | |
(loop | |
[i 5] | |
(if | |
(<= (sq i) n) | |
(if | |
(or (divides? i n) (divides? (+ 2 i) n)) | |
false | |
(recur (+ i 6))) | |
true)))) | |
(is (prime? 1020379)) | |
(is (not (prime? 1020378))) | |
(is (not (prime? 1)))) | |
(with-test | |
(defn | |
max-dividing-power | |
[p n] | |
(defn- | |
max-dividing-power-naive | |
[p n] | |
{:post [(natural? %)]} | |
(loop | |
[p-to-e 1 e 0] | |
(if (divides? p-to-e n) (recur (* p p-to-e) (inc e)) (dec e)))) | |
(defn- | |
find-power | |
[p-to-e e] | |
(+ e (max-dividing-power-naive p (quot n p-to-e)))) | |
(defn- | |
find-start | |
[p-to-e e] | |
(let | |
[p-to-e2 (sq p-to-e)] | |
(cond | |
(= p-to-e2 n) | |
(duple e) | |
(> p-to-e2 n) | |
(find-power p-to-e e) | |
(divides? p-to-e2 n) | |
(if | |
(divides? p (quot n p-to-e2)) | |
(find-start p-to-e2 (duple e)) | |
(duple e)) | |
:else | |
(find-power p-to-e e)))) | |
(cond | |
(one? p) | |
1 | |
(not (divides? p n)) | |
0 | |
:else | |
(assert (find-start p 1) natural?))) | |
(is | |
(= (max-dividing-power 3 27) 3) | |
(= (max-dividing-power 3 (duple 27)) 3))) | |
(with-test | |
(defn odd-prime? [n] (and (odd? n) (prime? n))) | |
(is (odd-prime? 3)) | |
(is (not (odd-prime? 2)))) | |
(with-test | |
(defn randint [a b] (+ (rand-int (+ b (- a) 1)) a)) | |
(is | |
(= | |
(vec | |
(sort | |
(keys | |
(frequencies | |
(take (pow 2 8) (repeatedly (fn* [] (randint 2 8)))))))) | |
(inrange 2 8)))) | |
(with-test | |
(defn- | |
prime-strong-pseudo-single? | |
"this thing maybe possibly works" | |
[n] | |
{:pre [(natural? n)]} | |
(cond | |
(<= 4 n) | |
(let | |
[a (randint 2 (- n 2)) g (gcd a n)] | |
(if | |
(< 1 g) | |
g | |
(loop | |
[ν 0 m (dec n)] | |
(cond | |
(even? m) | |
(recur (inc ν) (quot m 2)) | |
:else | |
(let | |
[b (modular-expt a m n)] | |
(cond | |
(one? b) | |
:probably-prime | |
:else | |
(loop | |
[i 0 b b b-old b] | |
(if | |
(and (< i ν) (not (one? b))) | |
(recur (inc i) (mod (sq b) n) b) | |
(if | |
(one? b) | |
(let | |
[g (gcd (inc b-old) n)] | |
(if (or (one? g) (= g n)) :probably-prime g)) | |
:composite))))))))) | |
(one? n 1) | |
:composite | |
:else | |
:probably-prime)) | |
(is (= (prime-strong-pseudo-single? 4) 2)) | |
(is (#{3 2} (prime-strong-pseudo-single? 6))) | |
(is (#{:composite 4 2} (prime-strong-pseudo-single? 8))) | |
(is (= (prime-strong-pseudo-single? 5) :probably-prime)) | |
(is (= (prime-strong-pseudo-single? 7) :probably-prime)) | |
(is (= (prime-strong-pseudo-single? 11) :probably-prime))) | |
(quote | |
(defn- | |
prime-strong-pseudo-explanation | |
[n] | |
(loop | |
[trials 1.0E7 result (prime-strong-pseudo-single? n)] | |
(cond | |
(zero? trials) | |
:very-probably-prime | |
(= result (quote probably-prime)) | |
(recur (dec trials) (prime-strong-pseudo-single? n)) | |
:else | |
result)))) | |
(quote | |
(defn | |
prime-strong-pseudo? | |
[n] | |
(let | |
[explanation (prime-strong-pseudo-explanation n)] | |
(or (= explanation :very-probably-prime) (= explanation true))))) | |
(quote | |
(defn | |
small-prime-factors-over | |
[n p] | |
{:pre [(natural? p)]} | |
(cond | |
(neg? p) | |
(raise-argument-error (quote small-prime-factors-over) "Natural" p) | |
(< n p) | |
(list) | |
(= n p) | |
(list (list p 1)) | |
(prime? n) | |
(list (list n 1)) | |
(divides? p n) | |
(let | |
[m (max-dividing-power p n)] | |
(cons | |
(list p m) | |
(small-prime-factors-over (quot n (expt p m)) (next-prime p)))) | |
:else | |
(recur n (next-prime p))))) | |
(quote | |
(with-test | |
(defn as-power [a]) | |
(are | |
[[a b e]] | |
(let [[b2 e2] (as-power a)] (and (= b b2) (= e e2))) | |
(apply (as-power a)) | |
[27 3 3] | |
[28 28 1] | |
[(* 5 5 7 7 7) (* 5 5 7 7 7) 1] | |
[(* 5 5 7 7 7 7) 245 2]))) | |
(quote | |
(with-test | |
(defn | |
perfect-power? | |
[a] | |
(and | |
(not (zero? a)) | |
(let [(base n) (as-power a)] (and (> n 1) (> a 1))))) | |
(are | |
[a b] | |
(= (perfect-power a) b) | |
3 | |
false | |
9 | |
true | |
(pow 12 7) | |
true | |
(dec (pow 12 7)) | |
false))) | |
(quote | |
(defn | |
powers-of | |
[a n] | |
(loop | |
[i 0 a-to-i 1] | |
(if (<= i n) (cons a-to-i (loop (+ i 1) (* a-to-i a))) (quote ()))))) | |
(quote | |
(with-test | |
(defn as-power [a]) | |
(are | |
[[a b e]] | |
(let [[b2 e2] (as-power a)] (and (= b b2) (= e e2))) | |
(apply (as-power a)) | |
[27 3 3] | |
[28 28 1] | |
[(* 5 5 7 7 7) (* 5 5 7 7 7) 1] | |
[(* 5 5 7 7 7 7) 245 2]))) | |
(quote | |
(with-test | |
(defn | |
perfect-power? | |
[a] | |
(and | |
(not (zero? a)) | |
(let [(base n) (as-power a)] (and (> n 1) (> a 1))))) | |
(are | |
[a b] | |
(= (perfect-power a) b) | |
3 | |
false | |
9 | |
true | |
(pow 12 7) | |
true | |
(dec (pow 12 7)) | |
false))) | |
(quote | |
(defn | |
powers-of | |
[a n] | |
(loop | |
[i 0 a-to-i 1] | |
(if (<= i n) (cons a (loop (+ i 1) (* a-to-i a))) (quote ()))))) | |
(successful? (run-tests)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment