## 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 xy 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) / (xn), 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 n2 = a b. That means each pair of divisors of n2 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 n2, 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.

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