Largest divisor such that two numbers divided by it round to the same value?

Question

I've got an algorithm that can be interpreted as dividing up the number line into an equal number of chunks. For simplicity, I'll stick with [0,1), it will be divided up like so:

0|----|----|----|----|1

What I need to do is take a range of numbers [j,k) and find the largest number of chunks, N, up to some maximum M, that will divide up the number line so that [j,k) still all fall into the same "bin". This is trickier than it sounds, as the range can straddle a bin like so:

    j|-|k
0|----|----|----|----|1

So that you may have to get to quite a low number before the range is entirely contained. What's more, as the number of bins goes up, the range may move in and out of a single bin, so that there's local minima.

The obvious answer is to start with M bins, and decrease the number until the range falls into a single bin. However, I'd like to know if there's a faster way than enumerating all possible divisions, as my maximum number can be reasonable large (80 million or so).

Is there a better algorithm for this?

Answer 1

Here I would like to give another heuristic, which is different from btilly's.

The task is to find integers m and n such that m / n <= j < k <= (m + 1) / n, with n as big as possible (but still under M).

Intuitively, it is preferable that the fraction m / n is close to j. This leads to the idea of using continued fractions.

The algorithm that I propose is quite simple:

1. calculate all the continued fractions of j using minus signs (so that the fractions are always approching j from above), until the denominator exceeds M;
2. for each such fraction m / n, find the biggest integer i >= 0 such that k <= (m * i + 1) / (n * i) and n * i <= M, and replace the fraction m / n with (m * i) / (n * i);
3. among all the fractions in 2, find the one with biggest denominator.

The algorithm is not symmetric in j and k. Hence there is a similar k-version, which in general should not give the same answer, so that you can choose the bigger one from the two results.

---

Example: Here I will take btilly's example: j = 0.6 and k = 0.65, but I will take M = 10.

I will first go through the j-procedure. To calculate the continued fraction expansion of j, we compute:

  0.6
= 0 + 0.6
= 0 + 1 / (2 - 0.3333)
= 0 + 1 / (2 - 1 / (3 - 0))

Since 0.6 is a rational number, the expansion terminates in fintely many steps. The corresponding fractions are:

0 = 0 / 1
0 + 1 / 2 = 1 / 2
0 + 1 / (2 - 1 / 3) = 3 / 5

Computing the corresponding i values in step 2, we replaces the three factions with:

0 / 1 = 0 / 1
1 / 2 = 3 / 6
3 / 5 = 6 / 10

The biggest denominator is given by 6 / 10.

---

Continue with the example above, the corresponding k-procedure goes as follows:

  0.65
= 1 - 0.35
= 1 - 1 / (3 - 0.1429)
= 1 - 1 / (3 - 1 / (7 - 0))

Hence the corresponding fractions:

1 = 1 / 1
1 - 1 / 3 = 2 / 3
1 - 1 / (3 - 1 / 7) = 13 / 20

Passing step 2, we get:

1 / 1 = 2 / 2
2 / 3 = 6 / 9
13 / 20 = 0 / 0 (this is because 20 is already bigger than M = 10)

The biggest denominator is given by 6 / 9.

---

EDIT: experimental results.

To my surprise, the algorithm works better than I thought.

I did the following experiment, with the bound M ignored (equivalently, one can take M big enough).

In every round, I generate a pair (j, k) of uniformly distributed random numbers in the inteval [0, 1) with j < k. If the difference k - j is smaller than 1e-4, I discard this pair, making this round ineffective. Otherwise I calculate the true result trueN using naive algorithm, and calculate the heuristic result heurN using my algorithm, and add them to statistic data. This goes for 1e6 rounds.

Here is the result:

effective round     =  999789
sum of trueN        =  14013312
sum of heurN        =  13907575
correct percentage  =  99.2262 %
average quotient    =  0.999415

The correct percentage is the percentage of effective rounds such that trueN is equal to heurN, and the average quotient is the average of the quotient heurN / trueN for all effective rounds.

Thus the method gives the correct answer in 99%+ cases.

I also did experiments with smaller M values, and the results are similar.

Answer 2

The best case for the bin size must be larger than k-j.

Consider the number line segments [0..j]and [k..1). If we can divide both of the partial segments into parts using the same bin size, we should be able to solve the problem.

So if we consider gcd((j-0)/(k-j), (1-k)/(k-j)), (where we use the greatest-integer-function after the division), we should be able to get a good estimate, or the best value. There are corner cases: if (k-j) > j or (k-j) > (1-k), the best value is 1 itself. So a very good estimate should be min(1, (k-j) * gcd((j-0)/(k-j), (1-k)/(k-j)))

Answer 3

Let's turn this around a bit.

You'd like to find m, n as large as you can (though n < M) with m/n close to but less than j and k <= (m+1)/n.

All promising candidates will be on the https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree. Indeed you'll get a reasonably good answer just walking the Stern-Brocot tree to find the last "large rational" fitting your limit just below j whose top is at k or above.

There is a complication. Usually the tree converges quickly. But sometimes the Stern-Brocot tree has long sequences with very small gaps. For example the sequence to get to 0.49999999 will include 1/3, 2/5, 3/7, 4/9, ... We always fall into those sequences when a/b < c/d, then we take the median (a+c)/(b+d) and then we walk towards one side, so (a+i*c)/(b+i*d). If you're clever, then rather than walking the whole sequence you can just do a binary search for the right power of i to use.

The trick to that cleverness is to view your traversal as:

1. Start with 2 "equal" fractions.
2. Take their median. If that exceeds M then I'm done. Else figure out which direction I am going from that.
3. Try powers of 2 in (a+i*c)/(b+i*d) until I know what range i is in for my range and M conditions.
4. Do binary search to find the last i that I can use.
5. (a+i*c)/(b+i*d) and (a+i*c+c)/(b+i*d+d) are my two new equal fractions. Go back to the first step.

The initial equal fractions are, of course, 0/1 and 1/1.

This will always find a decent answer in O(log(M)) operations. Unfortunately this reasonably good answer is NOT always correct. Consider the case where M = 3, j=0.6 and k=0.65. In this case the heuristic would stop at 1/2 while the actual best answer is 1/3.

Another way that it can fail is that it only finds reduced answers. In the above example if M was 4 then it still thinks that the best answer is 1/2 when it is actually 1/4. It is easy to handle this by testing whether a multiple of your final answer will work. (That step will improve your answer a fixed, but reasonably large, fraction of the time.)