Diophantine Reciprocals
September 19, 2014
We’ll take the second half of the question first: make a list of pairs that solve a particular n. It is obvious that both x and y must be greater than n, and we arbitrarily choose x ≤ y to prevent duplicates. Further, for any given n, the pair 2n,2n is always a solution, so we can find all solutions by iterating x from n + 1 to 2n. Solving for y gives y = (n * x) / (x − n), so we just look for integral y in a loop:
(define (xy-pairs n)
(let loop ((x (+ n 1)) (xys (list)))
(if (< (+ n n) x)
(reverse xys)
(let ((y (/ (* x n) (- x n))))
(if (integer? y)
(loop (+ x 1) (cons (list x y) xys))
(loop (+ x 1) xys))))))
Thus, to find the least n with 1000 solutions, say:
> (let loop ((n 1))
(if (<= 1000 (length (xy-pairs n))) n (loop (+ n 1))))
180180
That took about four hours on my slow machine at home.
A better solution to the counting problem uses a little bit of math. Since both x and y are greater than n, we can rewrite the original Diophantine equation as 1 / (n + a) + 1 / (n + b) = 1 / (n, then doing some algebra we get (n + b) / ((n + a) (n + b)) + (n + a) / ((n + a) (n + b)) = 1 / n, then n(2n + a + b) = ((n + a) (n + b)) and finally n 2 = a b. That means each pair of divisors of n 2 represents a solution.
From a previous exercise we have a function that computes the number of divisors of a number, which we can use to compute the count of solutions for a given n without enumerating the solutions, and we use a 2,3,5-wheel to perform the factorization:
(define (numdiv n)
(let ((fs (factors n)))
(let loop ((prev (car fs)) (fs (cdr fs)) (f 2) (d 1))
(cond ((null? fs) (* d f))
((= (car fs) prev) (loop prev (cdr fs) (+ f 1) d))
(else (loop (car fs) (cdr fs) 2 (* d f)))))))
(define (xy-count n) (/ (+ (numdiv (* n n)) 1) 2))
> (let loop ((n 1)) (if (<= 1000 (xy-count n)) n (loop (+ n 1))))
180180
That runs in about three-and-a-half minutes on my home computer, about two orders of magnitude faster.
Of course, if we want to find the factors of n 2, we can just find the factors of n and consider each one twice:
(define (numdiv2 n)
(let ((fs (factors n)))
(let loop ((prev (car fs)) (fs (cdr fs)) (f 2) (d 1))
(cond ((null? fs) (* d (+ f 1)))
((= (car fs) prev) (loop prev (cdr fs) (+ f 2) d))
(else (loop (car fs) (cdr fs) 2 (* d (+ f 1))))))))
(define (xy-count n) (/ (+ (numdiv2 n) 1) 2))
> (let loop ((n 1)) (if (<= 1000 (xy-count n)) n (loop (+ n 1))))
180180
That runs in about two seconds on that same slow machine, another two orders of magnitude improvement.
You can run the program at http://programmingpraxis.codepad.org/cvHMoejH. I would be very impressed by any candidate who stood at a white board and solved that problem in a one-hour interview. Melinda Lu gives a good description of the solution.
Programming Praxis 
Here’s another way to solve this (taking the problem to be finding the smallest number whose square has at least 2000 divisors).
A number of the form (p1^a1)..(pn^an), for primes p1,…pn, has (a1+1)…(an+1) divisors and its square has (2a1+1)…(2an+1) divisors. For a given a1…an, which we can assume to be in non-increasing order, the smallest such number is where p1 = 2, p2 = 3 etc. Therefore by generating all such smallest numbers in order and checking the number of divisors, we will find the solution we seek. So, we need to generate all non-increasing sequences a1,a2,a3,… in order of (2^a1)(3^a2), and we can do this by using a variant of Dijkstra’s shortest path algorithm on a graph of sequences, with edges (a1,…,an) -> (a1,…,an,1) and (a1,…,ai,…,an) -> (a1,…,ai+1,…an) (as long as that is non-increasing), so, for example, ()->(1)->(2)->(2 1)->(2 2)->(3 2)->(3 2 1) is a path in the graph, representing 1->2->4->12->36->72->360. We maintain a queue of unvisited sequences and each stage we output the smallest unvisited sequence/number and add its successors into the queue (if they aren’t already there).
Here’s a C++ implementation using STL heap operations to manage the queue and an STL set to keep track of seen items, runtime is about 4 ms.
#include <stdlib.h> #include <assert.h> #include <iostream> #include <vector> #include <set> #include <algorithm> const int K = 2; // Calculate divisors of n^K const int NPRIMES = 10; int primes[NPRIMES] = { 2,3,5,7,11,13,17,19,23,29 }; struct Entry { int n; // The represented number int ndivisors; // Number of divisors of interest std::vector<int> factors; // Coefficients of prime factors Entry() : n(1), ndivisors(1) {} // Comparison for heap operations static bool comp(const Entry *e1, const Entry *e2) { return e1->n > e2->n; } // Attempt to make a successor node by incrementing element k // Return NULL if invalid or already seen const Entry *extend(int k, const std::set<int> &seen) const { assert(k < NPRIMES); int n1 = n*primes[k]; // Check we haven't seen the node before. if (seen.count(n1) > 0) return NULL; int fsize = factors.size(); // Check result is non-increasing sequence if (k > 0 && k < fsize && factors[k-1] == factors[k]) return NULL; // All's well, make the new entry. Entry *e = new Entry(); e->n = n1; e->factors = factors; if (k == fsize) e->factors.push_back(0); int t = e->factors[k]; e->ndivisors = ndivisors/(1+K*t)*(1+K+K*t); e->factors[k] = t+1; return e; } }; int main(int argc, char *argv[]) { int N = atoi(argv[1]); std::set<int> seen; std::vector<const Entry*> q; q.push_back(new Entry()); while(true) { std::pop_heap(q.begin(),q.end(),Entry::comp); const Entry *e = q.back(); q.pop_back(); if (e->ndivisors > N) { std::cout << e->n << "\n"; break; } for (int i = 0; i <= (int)e->factors.size(); i++) { const Entry *e1 = e->extend(i,seen); if (e1 != NULL) { q.push_back(e1); seen.insert(e1->n); std::push_heap(q.begin(),q.end(),Entry::comp); } } } }Erlang solution, factoring is done by a 2,3,5,7 factor wheel.
-mode(compile). -import(lists, [foldl/3]). divisors_of_sq(N) -> foldl(fun ({_, Exp}, A) -> (Exp*2 + 1) * A end, 1, factor:factorize(N)). solution_count(N) -> (divisors_of_sq(N) + 1) bsr 1. solve(Limit) -> solve(1, Limit). solve(Start, Limit) -> N = solution_count(Start), if N > Limit -> Start; true -> solve(Start + 1, Limit) end. main(_) -> Euler108 = solve(1000), io:format("1/x + 1/y = 1/~w has ~w solutions.~n", [Euler108, solution_count(Euler108)]).Here’s another solution: if n = (p^k)m where p is prime and doesn’t divide m, then divisors(n) = (k+1)*divisors(m), with m < n, so for each n, we just need to apply trial division until we have found a single prime power factor, then apply the recurrence with a table lookup:
#include
#include
const int K = 2; // We want divisors of n^2
int divisors(int n, const std::vector &a)
{
// Find a prime power factor, reduce n and
// apply recurrence.
int k = 0;
for (int p = 2; p*p 0) return (1+K*k)*a[n];
}
return 1+K; // n is prime
}
int main()
{
int N = 2000;
std::vector a;
a.push_back(0); a.push_back(1);
for (int n = 2; ; n++) {
a.push_back(divisors(n,a));
if (a[n] >= N) {
std::cout << n << "\n";
break;
}
}
}
Don’t know what happened to the formatting there, let’s try again:
#include <vector> #include <iostream> const int K = 2; // We want divisors of n^2 int divisors(int n, const std::vector<int> &a) { // Find a prime power factor, reduce n and // apply recurrence. int k = 0; for (int p = 2; p*p <= n; p++) { while (n%p == 0) { k++; n/=p; } if (k > 0) return (1+K*k)*a[n]; } return 1+K; // n is prime } int main() { int N = 2000; std::vector<int> a; a.push_back(0); a.push_back(1); for (int n = 2; ; n++) { a.push_back(divisors(n,a)); if (a[n] >= N) { std::cout << n << "\n"; break; } } }