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[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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s

%d bloggers like this: