Nearly-Square Divisors
August 5, 2016
We have today a fun little problem from number theory:
Given positive integers n, a and b such that n = a · b with a ≥ b, find a and b such that the difference a − b is as small as possible.
For n square, the solution is just the square root of n; for instance, with n = 36, a = b = 6. Otherwise, a and b will be the two divisors nearest the square root; for instance, with n = 60, a = 10 and b = 6.
Your task is to write a program to find a and b as described above; use your program to find the nearly square divisors of n = 224403121196654400. 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.
Programming Praxis 
Another method is to start with the square root of n, rounded down, then decrement by one until a divisor is found. On my computer this took 0.6 sec.
The pitfall with any method is to deal with a prime n.
In Python. The divisors are scanned for the maximum value less equal isqrt(n).
def nsd(n): s = isqrt(n) a = max(d for d in divisors(n) if d <= s) return a, n // a@Ernie had the same idea! Here’s a solution that does that in Racket.
(define (near-squares n) (let loop ((a (floor (sqrt n)))) (let ((b (/ n a))) (if (and (integer? b) (= n (* a b))) (list (inexact->exact a) (inexact->exact b)) (loop (sub1 a))))))I came up with the same idea as the previous commenters. For large n with a higher number of divisors this may become incredibly time consuming, e.g. n = 100#. So I guess there’s a quicker solution than simply testing all possible numbers from floor(sqrt(n)) downward.
Python:
import math def nsd(n): s = (int)(math.sqrt(n)) while n % s: s -= 1 return s, n // s print (nsd(224403121196654400))Output in approx. 0.012s: (473670855, 473753280)
Anyhow, this solution is based on the slow brute-force method. I’m almost certain there’s a much faster number-theoretical way to solve this for huge n.
@Ernie and @Milbrae. The number N=224403121196654400 has an incredibly large number of divisors and is the worst-case scenario for a solution using divisors and the best-case scenario for the solution walking down from the root. Now try to use both methods in the range [N-10, N+10). Timing is 90 ms for the divisor solution and the other solution takes forever.
Paul: That is point I was trying to make in my first post: if N has very few divisors, decrementing from the sqrt may take forever, and if N is prime, well don’t ask.
224403121196654400 is the 22nd “colossally abundant” number (and mentioned here: https://programmingpraxis.com/2016/04/22/gcd-sum/#comment-59608)
Finding nearly-square divisors seems about as hard as factoring (if we have a method for finding them we can apply it recursively to get a full factorization, and if we have a factorization, finding such divisors is straightforward).
Fermat factorization seems relevant too, but maybe only applies to odd numbers.
@Paul: You may have a point there. Anyhow, consider n having at least 30 distinct prime factors, each appearing once. n will then have 2**30 factors. You wouldn’t want to try almost all of them to find the correct NSD.
There’s a similar problem on Project Euler (https://projecteuler.net/problem=266). I don’t think the people you have already solved this did it through simply factoring.
@Milbrae: Good point, actually it’s a variant on the knapsack problem – find a maximal subset (or sub-multiset) of prime factors of n, the sum of whose logs is less than log n / 2 (and there are ways of solving the knapsack problem that are better than testing each alternative).
@Milbrae. Thanks for pointing me to Euler 266. I solved for the product of first 30 primes with divisors and it took about 4 hours. The Euler problem needs the first 42 primes. That will indeed need another approach.
@Matthew: A varint of knapsack? Could you elaborate on this?
@Milbrae: One form of knapsack problem is: given a set of real numbers S, and a limit k, find the maximal subset sum less (or equal) to k (ie. find a subset with sum less than or equal to k, and greater or equal to the sum of any other such subset). This is exactly what we want for the nearly-square root of a product of single primes (as in the Project Euler problem) – S is the set of logs of the primes, and k is the log of the square root.
Here’s a simple Haskell program that computes the log of the nearly square root of the product of the first 20 primes:
-- k target s: return largest subset sum of s <= target k target [] = 0 k target (a:s) = if a > target then k target s else max (a + k (target-a) s) (k target s) primes = [2,3,5,7,11,13,17,19,23,29, 31,37,41,43,47,53,59,61,67,71] s = map log(reverse primes) target = (sum s)/2 r = k target s main = print target >> print (k target s)Scan the tree of subsets, discarding a branch which cannot lead to a solution (here just subsets that would be too large – we could also keep track of the best solution found and discard branches that cannot lead to a better solution).
@Matthew: Wow!! Works fine and extremely fast. I assume correctly the second line of the output is the exponent to e to get the result? I know practically nothing about Haskell.
[…] a recent exercise, we discussed the problem of computing the nearly square divisors of a composite number n = a · b, […]
[…] nearly square divisors exercise has generated a considerable amount of interest, and several excellent solutions in the comments. […]
package test;
import java.util.*;
public class SquareDivisors {
public static int squareDivisor(int n){
int diff=10000;
if (n < 0){
return diff;
}
List pairs = new ArrayList();
//no need to look beyond sqrt hence check for (i*i )<=n
for(int i=1;(i*i)<=n;i++){
if (n%i==0){
System.out.println("Pairs : "+i +" "+ (n/i));
pairs.add(new Pairs(i,n/i));
}
}
for(Pairs pair: pairs){
diff= (diff < (Math.abs(pair.a-pair.b)))?diff:(Math.abs(pair.a-pair.b));
}
return diff;
}
public static void main(String args[]){
System.out.println("Diff : "+squareDivisor(36));
}
}
class Pairs{
public Pairs(int a,int b){
this.a=a;
this.b=b;
}
public int a,b;
}