How to find a pair of divisors as close as possible to each other?

Question

For a given integer $n\in\mathbb{N}^*$, I want to find a pair $(x,y)\in{\mathbb{N}^*}^2$ such that $x*y=n$ and $|y-x|$ is as small as possible.

A naive algorithm I found is :

for k = floor(sqrt(n)) to 1 step -1
  if n/k is integer then
    return (k, n/k)
  end if
end for

Although I am almost certain there is a more efficient way to do this, I cannot find one.

Answer 1

Note that although there could be more efficient algorithms, what you found is not a bad algorithm. I will analyze this a little bit.

Let us write down the problem formally as a constrained discrete optimization problem: $$ (x^\star, y^\star) = \arg\min_{(x,y) \in \mathbb{Z_{+}}} \Big(E \triangleq \big(y-x\big)^2 + \lambda(x*y - n)^2 \Big) $$ Note that for positive integers the $L_2$ norm will attain the same minimum as $L_1$. So I'll stick to it. $\lambda$ is a Lagrange multiplier and the right-most term is the non-linear constraint. Even if we can compute the partial derivatives and use the method of Lagrange multipliers, the problem is not easy to solve due to the parameters restricted to a discrete space. Hence we will relax the integer constraint to real numbers and solve the problem for $\mathbb{R}$. The relaxed optimum would be attained as follows: $$ x^\star = \arg\min_{x \in \mathbb{R}} \Big(\tilde{E} \triangleq \big(x-\frac{n}{x}\big)^2 \Big) $$ We achieved this by substitution the constraint : $y \gets \frac{n}{x}$. Note that this is also suggested by the Largange optimization: $$ \frac{\partial E}{\partial \lambda} = xy - n^2 = 0 \implies x=\frac{n}{y} \text{ or } y=\frac{n}{x}. $$

Now differentiating w.r.t. $x$: $$ \begin{align} \frac{d \tilde{E}}{dx} = 2\big(x-\frac{n}{x}\big) &= 0\\ x^2 &= n \\ x = \sqrt{n} \end{align} $$ gives the optimal solution for the relaxed problem. For negative numbers $x = -\sqrt{n}$ is also the optimum. This we could already detect from the sign of the input. Having solved the relaxed problem, a typical optimizer would operate by projecting the solution into the feasible set and/or applying a projected gradient descent. The initial feasible solution is $x_0 = \lfloor \sqrt{n} \rceil$. $\lfloor \cdot \rceil$ is the rounding operation or in fancy words the projection onto the integers. From here, it is possible to walk on the Euclidean discrete manifold of 1-dimensional integers by your method: taking steps of $+1$ or $-1$ depending on the direction. That would amount to a gradient descent by the step size $h=1$.

And these steps are what you are following implicitly when arriving at the algorithm you proposed. Unfortunately, the worst case complexity of this algorithm is attained when the number is prime and we have only 2 divisors: $1$ and $n$, the number itself. Therefore, the worst case complexity reaches to $O(\sqrt{n})$. As $n\rightarrow \infty$, this would resemble a linear search (it is worse than $\log$). Now, the set of integers is a totally ordered set. Using this property, instead of using a linear search, one could think of a binary search variant. But as clarified by @WolfgangBangerth, such an algorithm would not exist bounding the complexity in: $O(\log(\sqrt{n}))<C\leq O(\sqrt{n})$.

Even though this might not be the exact answer you were looking for, I thought it would give a perspective to develop better solutions. For instance, the integer programming community might have a better solution.

Answer 2

As @FedericoPoloni commented on the post above, a natural way to expedite this search (relative to simple descending trial division) is by Fermat's Factorization Algorithm. As the Wikipedia article illustrates, the Fermat approach can take bigger steps than trial division, but at a cost of more arithmetic (square roots in extended precision versus only integer division).

Fermat's idea is to search for a difference of integers squares $N = a^2 - b^2$, so that $N = (a + b)(a - b)$. If the search starts with $a \ge \sqrt N$, then incrementally increasing $a$ finds the first such pair of "cofactors" of $N$, $c = a+b, d = a-b$ as close together as possible. CAVEAT: The cofactor pair $c,d$ would then have the same parity!

This means that we may either (i) not be able to express $N$ as a difference of two integer squares (in the case that $N = 4k+2$) or (ii) find a pair of same parity cofactors which are not the desired closest together pair of cofactors. An example of the first problem is $N=6$, which is not a difference of integer squares as shown in answers to this older Math.SE Question.

An example of the second problem is $N=12$, which can be expressed as a difference of squares $12 = 4^2 - 2^2 = 6\cdot 2$, but whose closest cofactors are instead $12 = 4\cdot 3$ (and not of the same parity).

Of course if $N$ is odd, it has only odd divisors, so any two cofactors will be the same parity (odd). To handle the case when $N$ is even, one can apply Fermat's method to $4N$. Now if we locate a cofactor pair of the same parity, they will both be even! So we will have $4N = (2c)(2d)$, and hence $N = cd$. Furthermore our search will locate the closest same parity cofactors $2c,2d$ of $4N$. It follows that $c,d$ are the closest cofactor pair of $N$ (though not necessarily of the same parity).

If $N$ is odd, then while a correct result will obtain from applying Fermat to $4N$, there is an advantage to applying directly to $N$. The oddness of $N$ means that numbers of the form $4k+1$ will be expressed as an odd square minus an even square, so our search is limited to checking $a$ odd. On the other hand numbers of the form $4k-1$ will be an even square minus an odd square, and again our search is limited (to checking $a$ even).

Algorithm: Fermat's Factorization Method (FFM)

Given as input positive integer $N \ge 1$ which is odd or a multiple of $4$:

0. Initialize integer $a = \lceil \sqrt N \rceil$. If $a = \sqrt N$, stop. OUTPUT: The closest cofactors are $c = d = \sqrt N$. Otherwise:

1. Now $a^2 > N$, so compute $a^2 - N$ and $b = \lfloor \sqrt{a^2 - N} \rfloor$. If $b = \sqrt{a^2 - N}$, stop. OUTPUT: The closest (same parity) cofactors are $c = a + b, d = a - b$. Otherwise:

2. Increment $a$ suitably (parity consideration) and go to Step 1.

---

It is possible to analyze at what point trial division becomes as efficient as Fermat's method, due to its trade-off of step size versus complexity. I plan to do an extended precision implementation, with the speculation that the time required for Fermat's step is about three times the time needed for trial division.