## 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 2*n*,2*n* is always a solution, so we can find all solutions by iterating *x* from *n* + 1 to 2*n*. 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*(2*n* + *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.

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.

Erlang solution, factoring is done by a 2,3,5,7 factor wheel.

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: