Diophantine Reciprocals
September 19, 2014
Career Cup claims that Amazon asked this as an interview question; it is also Problem 108 at Project Euler:
In the following equation x, y and n are positive integers: 1 / x + 1 y = 1 / n. For n = 4 there are exactly three distinct solutions: 1/5 + 1/20 = 1/6 + 1/12 = 1/8 + 1/8 = 1/4. What is the least value of n for which the number of distinct solutions exceeds one thousand?
Your task is to solve Amazon’s question; you might also like to make a list of the x, y pairs that sum to a given n. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
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; } } }