## 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.

Pages: 1 2

### 4 Responses to “Diophantine Reciprocals”

1. matthew said

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);
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);
}
}
}
}
```
2. David said

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)]).
```
3. matthew said

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;
}
}
}

4. matthew said

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;
}
}
}
```