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.

Advertisements

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: