I’ve had to hit the pause button on the blog for a bit because life got in the way, but I haven’t stopped my mathematical explorations. Actually, I got turned on to Project Euler in the meantime, which became something of an addiction. Most people solve the problems in C/C++ or Python, but I decided to take the opportunity to sharpen my skills in Clojure (which is, let’s face it, the best-designed programming language out there 😉).
In the beginning, the “mathematical discovery factor” increased with the difficulty, but after solving about 90 problems, those returns have plateaued and are now declining in relation to the diffficulty. So now it’s time to start posting about some of the cool things I’ve discovered along the way, and use this blog to explore those topics further.
Prime time
Thanks to Chapter 1 of A Programmer’s Introduction to Mathematics, I had already developed an interest in number theory and cryptography. This interest was further ignited by Project Euler, in whose problems number theory factors heavily (pun intended).
Number theory is basically the study of the integers, \(ℤ\). This in turn shines the spotlight on prime numbers, a very special category of the integers due to their many interesting properties.
In order to do anything interesting in number theory, we need to find a way to test the primality of a number and to generate a sequence of prime numbers.
In imperative languages, it’s typical to separate these two tasks, and approach the latter by implementing what is called a “prime sieve,” such as the Sieve of Eratosthenes, which walks up an infinite sequence of numbers and crosses off more and more composite numbers with each iteration. You might start with an enormous array of Trues and switch the indices of composite numbers to False as you go along.
Clojure, on the other hand, is a functional language, in which mutating sequences in-place is neither idiomatic nor clean. Lazy (lazily evaluated) infinite sequences are the way to go.
As the above link shows, it is possible to implement a Sieve of Eratosthenes in Clojure, but it’s quite difficult to read, and in any less than expert hands, quite ugly as well. Meanwhile, I’ve found that it’s simpler to filter primes in Clojure than sieving them, which also conveniently integrates the task of testing primality as well.
This approach is highly readable and, when memoized, its performance does not catastrophically degrade until you start needing primes greater than 1 million or so. Obviously, actual cryptographic applications would require a more industrial-strength implementation, as real-world cryptography deals with integers that could be as large as 256 bits (hundreds of digits long), but for armchair explorations of number theory, this approach will suffice.
Primes on trial
This method is called “trial division,” since it involves dividing \(n\) by numbers in a range to test for divisibility. First, we have to define divisibility. A dividend \(n\) is divisible by a divisor \(d\) if no remainder is left after the division, or in other words, if \(n\) is a multiple of \(d\).
In terms of modular arithmetic, this is conveniently expressed as a congruence: \(n \equiv 0 \mod d\).
(defn divisible? [n d]
(zero? (mod n d)))
The most naïve, brute-force approach is to literally divide \(n\) by all numbers in \([1, n)\) and return true (\(n\) is prime) if \(n\) is divisible by none of those numbers:
(defn brute-force-prime? [n]
(not-any? (partial divisible? n) (range 2 n)))
This is extremely slow, though. For a number like \(45565962173\), this is completely intractable, but even generating primes up to \(100000\) is impractical with this method.
(time (doall (filter brute-force-prime? (range 1 100000))))
;; Elapsed time: 67653.900002 msecs
A smarter way to do this starts with noticing that only potential divisors up to \(\sqrt{n}\) need to be tried, because any divisor less than that will have a complementary divisor on the other side of \(\sqrt{n}\). For example, \(\sqrt{28} \approx 5.3\). \(28\) is divisible by \(2\), and this division yields the complementary divisor \(14\). Likewise, \(4\) yields \(7\).
Integer square root
An upper integer bound \(\lfloor \sqrt{n} \rfloor\) capable of handling fairly large numbers can be implemented quite simply as (bigint (Math/sqrt n)), taking advantage of Java interop.
This has been effective for all of the Project Euler problems I’ve solved so far, though it is true that a proper arithmetical solution would not involve floating-point numbers. For sufficiently large \(n\), rounding errors would yield an upper bound that is too low to find all possible factors, especially in the worst case (where a very large \(n\) is the product of two very large primes that are similar but not equal in value).
I use the clojure.math.numeric-tower library here as
using Java interop for exponentiation ((Math/pow n e)) similarly introduces floating-point numbers, and (tower/expt n e) is ever-so-slightly faster than (reduce *' (repeat e n)).
Here is a bulletproof integer square root function, which uses the simpler, faster method for \(n < 10^24\) and a more sophisticated algorithm (which I don’t really understand) for larger numbers.
(require '[clojure.math.numeric-tower :as tower])
(defn isqrt
"floor(√n). When incremented, provides an upper bound for factorization."
;; Java interop is super fast but not accurate for n > 1E24 (approx) due to
;; floating-point rounding. Uses a slightly slower but pinpoint-precise method for n > 1E24.
[n]
(if (< n 1E24)
(-> (Math/sqrt n) bigint)
;; https://cs.stackexchange.com/a/30383
(let [half-bit-length (quot (.bitLength (bigint n)) 2)]
(loop [a (tower/expt 2 half-bit-length)
b a
c (*' a a)]
(cond
(zero? b) a
(> c n) (recur (-' a b) (quot b 2) (+ c (*' -2 a b) (*' b b)))
:else (recur (+' a b) (quot b 2) (+ c (*' 2 a b) (*' b b))))))))
Testing primality
Next, let’s implement some basic checks to eliminate the need to do any calculations in the majority of cases. When testing for primality, we generally only consider \(ℤ\), the positive integers. Furthermore, \(1\) is not considered prime. So \(n \leq 1\) is false right off the bat. Next, all primes are odd except \(2\), so make an exception for \(2\) and return false for all even? numbers.
Clojure’s built-in even? function has the added bonus of throwing an exception for Ratios and floats.
This leaves all odd numbers \(\geq 3\). Weeding out all even numbers means that there is no need to check for even divisors (odd numbers only have odd divisors).
The most idiomatic way to write this would be:
(cond
(<= n 1) false
(= n 2) true
(even? n) false ;; Will also weed out non-integers
:else (not-any? #(divisible? n %) (range 3 (inc (isqrt n)) 2)))
but this is actually about 1.5x as slow as a more verbose translation using loop/recur:
(defn naive-prime? [n]
(cond
(<= n 1) false
(= n 2) true
(even? n) false ;; Will also weed out non-integers
:else (let [lim (int-root n)]
(loop [i 3]
(cond
(>= i lim) true
(divisible? n i) false
:else (recur (+ i 2)))))))
(def prime? (memoize naive-prime?))
whose performance, non-memoized and memoized, is leaps and bounds above the original brute-force trial division function:
(require '[criterium.core :as c])
(c/quick-bench (naive-prime? 45565962173))
;; Execution time mean : 16.516873 ms
;; Execution time std-deviation : 1.053979 ms
(c/quick-bench (prime? 45565962173))
;; Execution time mean : 340.284934 ns
;; Execution time std-deviation : 10.658400 ns
Leveraging the primality test
This may not exactly be an industrial-strength method, but for my current purposes, it arguably allows more idiomatic and readable ways to accomplish tasks such as filter, take, and nth.
(c/quick-bench (doall (filter naive-prime? (range 1 100000))))
;; Execution time mean : 265.236965 ms
;; Execution time std-deviation : 6.784358 ms
(c/quick-bench (doall (filter prime? (range 1 100000))))
;; Execution time mean : 83.140024 ms
;; Execution time std-deviation : 1.655784 ms
(c/quick-bench (take 500 (filter prime? (range))))
;; Execution time mean : 88.645751 ns
;; Execution time std-deviation : 27.586353 ns
;; (dec n) gets the nth prime, because the sequence is zero-indexed
(c/quick-bench (nth (filter prime? (range)) (dec 500)))
;; Execution time mean : 2.473162 ms
;; Execution time std-deviation : 123.136771 µs
Prime factorization
With that out of the way, we come to the cornerstone of number theory:
All integers greater than 1 are either a prime number or can be expressed as a unique product of prime numbers.
The algorithm to find the prime factors of an integer \(n\) is surprisingly simple to implement. Start with an empty list of factors and an infinite list of primes. Begin with the first prime (that is, \(2\)) and divide \(n\) by \(2\) until you can’t anymore. Each time you divide, add a \(2\) to your list of factors.
Then, proceed to the next prime (\(3\)) and repeat up the list, adding factors, until the result of your division is \(1\).
Let’s try \(168 = 2 \times 2 \times 2 \times 3 \times 7\).
(loop [n 168
primes (filter prime? (range))
factors []]
(let [k (first primes)]
(cond
(= 1 n) factors
(divisible? n k) (recur (/ n k) primes (conj factors k))
:else (recur n (rest primes) factors))))
;; [2 2 2 3 7]
Checks out. You can try a bunch of random numbers and you’ll instantly get its prime factorization, which is something like a fingerprint, since it’s unique for every number. Pretty cool!
But again, a ridiculously large number like \(245454537724879\) is too much for a simple algorithm. Actually, it’s not so much that the number is too large, but that it has very large prime factors, so it takes a very long time to iterate that far up the list of primes.
This can be worked around somewhat if we find the factors of \(n\) first, then filter the prime? ones.
Regular factorization
This will look a bit like our initial foray into primality testing. Rather than checking if \(n\) has no (that is, not-any?) divisors \(1 < d < n\), we just filter all divisors \(1 \leq d < n\).
(defn brute-force-factors [n]
(filter (partial divisible? n) (range 1 n)))
As you can imagine, since this is running through every possible number, this will take way too long for large \(n\). We can use the same upper bound as before, although that means that for every positive result, we have to add not only \(d\), but \(\frac{n}{d}\) as well.
Here is a clean way of doing that:
(defn factors [n]
(->> (range 1 (inc (isqrt n)))
(filter #(divisible? n %))
(mapcat (fn [d] [d (/ n d)]))
(into (sorted-set))))
More optimized prime factorization
With this, we can now improve the performance of our prime factorization function by checking the size of the number first. If \(n\) is smaller than some threshold (let’s say 1 million), then iterate up an infinite list of primes as before. If \(n\) is larger than the threshold, though, then filter the prime? factors first, and iterate up that list instead.
(defn prime-factorization [n]
"Returns the prime factorization of an integer, e.g., 168 -> [2 2 2 3 7]."
(when (and (integer? n) (pos? n)) ;; Sanity check
(loop [n n
primes (if (> n 1000000)
(filter prime? (factors n))
(filter prime? (range)))
factors []]
(let [k (first primes)]
(cond
(= 1 n) factors
(divisible? n k) (recur (/ n k) primes (conj factors k))
:else (recur n (rest primes) factors))))))
This isn’t a perfect heuristic, as large numbers do not necessarily have large prime factors. For example, 12 bazillion-gajillion (\(12 \times 10^{??}\)) still has only \((2, 3, 5)\) as prime factors, just like \(120\). But at least it can reduce the time required for difficult cases, such as \(n = 23897538974893789\), and make them tractable.
(time (prime-factorization 23897538974893789))
;; Elapsed time: 30579.928312 msecs
;; [211 23357 4849016507]
(The brute-force version didn’t finish even after 5 minutes, so I abandoned it.)
Prime omega functions
There are a couple arithmetic that look and sound really impressive but are trivial to implement in code once you can factor an integer into primes.
One is the little omega function \(\omega(n)\), which counts the number of distinct prime factors of \(n\).
(count (distinct (prime-factorization n)))
The other is the big omega function \(\Omega(n)\), which counts the number of prime factors of \(n\) with multiplicity (i.e., if \(2\) appears more than once, count both instances).
(count (prime-factorization n))
Prime power representation
Given the prime factorization of a number, we can use this information to formulate its unique representation as a product of prime powers.
$$ n = p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}
= \prod_{i=1}^{k}p_{i}^{e_{i}} $$
Returning to our simple example \(168 = 2 \times 2 \times 2 \times 3 \times 7\) above, this can be written more succinctly as
$$ 168 = 2^3 \times 3 \times 7 $$
which is also called its canonical representation or standard form.
If we fill in the missing primes above with \(p_{i}^{0}\), which doesn’t affect the final product, then the sequence of \(e_1 \cdots e_k\) in the above notation can also be extracted from the prime factorization, which gives us an exponent vector (or unique prime signature). For \(168\), this would be
$$ 168 = 2^{\color{red}{3}} \times 3^{\color{red}{1}} \times 5^{\color{red}{0}} \times 7^{\color{red}{1}}
= \begin{bmatrix}3 & 1 & 0 & 1\end{bmatrix}$$
Here’s a simple way to implement that:
(defn prime-powers-integer [n]
"Returns the exponent vector for the prime power representation of an integer,
e.g., 168 = 2*2*2*3*7 = 2^3 * 3^1 * 5^0 * 7^1 -> (3 1 0 1)"
(when (and (integer? n) (pos? n))
(if (= n 1) [0]
(let [pf (prime-factorization n)
primes (filter prime? (range (inc (peek pf))))
freqs (frequencies pf)
get-exp (fn [p] (if-let [exp (get freqs p)]
exp 0))]
(map get-exp primes)))))
For \(n = 168\), the procedure works like this:
pf => [2 2 2 3 7]
primes => [2 3 5 7] ;; one of each prime up to the last prime in pf
freqs => {2 3, 3 1, 7 1}
(map get-exp primes) => (map get-exp [2 3 5 7]) => (3 1 0 1)
This works fine when the prime factors are small, but it isn’t the cleanest approach. Namely, one list of primes may already be generated by (prime-factorization n), so generating another (for the purpose of determining that \(87\) is the \(i\)th prime, for example) is not very efficient. If we limited the prime-factorization function to the iterative approach, we could generate a map containing {:i 2 :p 3 :e 3} for each prime, which is a starting point for generating both the canonical representation and the exponent vector.
However, since this is mostly for curiosity’s sake than actual applications, I’m leaving it at that. For cryptographic applications, what is more important than generating exponent vectors are:
- Finding the factors (prime or otherwise) of an integer themselves, and
- Specialized techniques for factoring specific types of numbers, such as semiprimes.
Now, for completeness’ sake, let’s write a function to convert an exponent vector back to an integer.
(defn prime-powers->num [pp]
(let [primes (filter prime? (range))]
(reduce *' (map tower/expt primes pp))))
(prime-powers->num [3 1 0 1])
;; 168
Negative exponents too?
While writing this, I learned that allowing negative exponents in prime factorizations enables you to represent not only all the integers (\(n \in ℤ\)), but all the rationals (\(n \in ℚ\)) as well.
For some rational number \(q = \frac{a}{b}\), let \(\vec a, \vec b\) be the exponent vectors for the integers \(a, b\). The exponent vector for \(q\) is then \(\vec q = \vec a - \vec b\). For example,
$$\begin{gathered}
q = \frac{171}{98} \\[0.8em]
\begin{aligned}
a = 171 &= 3^2 \times 19^1 \\
b = 98 &= 2^1 \times 7^2 \\[0.8em]
\end{aligned} \\
\begin{aligned}
\vec a &= \begin{bmatrix}\; \; \; 0 & 2 & 0 & \; \; \; 0 & 0 & 0 & 0 & 1\end{bmatrix} \\
\vec b &= \begin{bmatrix}\; \; \; 1 & 0 & 0 & \; \; \; 2 & \textcolor{#bbb}{0} & \textcolor{#bbb}{0} & \textcolor{#bbb}{0} & \textcolor{#bbb}{0}\end{bmatrix} \\
\vec q &= \begin{bmatrix}-1 & 2 & 0 & -2 & 0 & 0 & 0 & 1\end{bmatrix} \\[0.8em]
\end{aligned} \\
\begin{aligned}
q &= 2^{-1} \times 3^2 \times 7^{-2} \times 19^1 \\
&= \frac{1}{2} \times 9 \times \frac{1}{49} \times 19
\end{aligned}
\end{gathered}$$
The following code is something I scratched up quickly as Clojure does not have a simple way to implement something like map-longest. That is, a way to map over multiple collections, filling in dummy values to make each collection the same size as the largest—in this case, the light gray \(\textcolor{#bbb}{0}\) in \(\vec b\) above.
(defn prime-powers-rational [q]
(let [r (rationalize q)]
;; Sanity check to accept decimal representations of rational numbers
;; while still rejecting irrational numbers
(when (or (ratio? q) (= (double r) (double q)))
(let [n (prime-powers-integer (numerator r))
d (prime-powers-integer (denominator r))
pn (concat n (repeat (count d) 0))
pd (concat d (repeat (count n) 0))]
(->> (map - pn pd)
reverse
(drop-while zero?) ;; truncate zeros from the end
reverse)))))
(prime-powers-rational 171/98)
;; (-1 2 0 -2 0 0 0 1)
(prime-powers->num [-1 2 0 -2 0 0 0 1])
;; 171/98
Pretty cool!
While the rational version doesn’t reside strictly within the confines of number theory, as it goes beyond the integers, it is still mathematically interesting and can be wrapped up neatly in a bulletproof prime-powers function that can handle both integer and rational inputs.
(defn prime-powers [n]
(when (pos? n)
(if (integer? n)
(prime-powers-integer n)
(prime-powers-rational n))))
Apparently, this can be extended to rational powers of rational numbers as well, which would seem to allow irrational numbers to be represented as well, though that is starting to get a little deep for me 😅
That’s all for now. Stay tuned for more posts in this series!
References
- Fundamental theorem of arithmetic, Wikipedia
- Prime factorization, OEIS Wiki