Nearly Square Divisors, Revisited

August 23, 2016

In a recent exercise, we discussed the problem of computing the nearly square divisors of a composite number n = a · b, with a ≥ b, such that the difference ab is as small as possible. We gave two solutions: The first solution enumerated the divisors one-by-one, by trial division counting from 1, until reaching the square root, which is fast if n is small but slow in general. The second solution factored n, computed the divisors, the picked the two divisors in the middle of the list, which is fast in general but slow if n is highly composite and thus has a large number of divisors.

In a comment on that exercise, Matthew Arcus, who is a regular reader and commentor at Programming Praxis, gave a splendid answer to the problem; it relies on factoring n, but is fast even if n has a large number of divisors. His algorithm reduces the multiplication of the factors to addition of their logarithms, which means that the knapsack algorithm can be used to find the greatest sum less than the logarithm of the square root.

Your task is to implement Matthew’s algorithm to find the nearly-square divisors of a number. 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

6 Responses to “Nearly Square Divisors, Revisited”

  1. Tadeusz said

    How do this work? With logarithm?

  2. Gambiteer said

    You can use the knapsack algorithm with products, as well as with sums.

    On a 64-bit machine all numbers are fixnums, so things go even a bit faster:

    > (time (nsd-new 224403121196654400))
    (time (nsd-new 224403121196654400))
    63 ms real time
    64 ms cpu time (64 user, 0 system)
    no collections
    3280 bytes allocated
    no minor faults
    no major faults
    (473753280 473670855)
    > (time (nsd 224403121196654400))
    (time (nsd 224403121196654400))
    148 ms real time
    148 ms cpu time (148 user, 0 system)
    77 collections accounting for 62 ms real time (56 user, 0 system)
    165921824 bytes allocated
    no minor faults
    no major faults
    (473753280 473670855)

  3. Gambiteer said

    And here’s some code (without sum or factors):

    (define (multiply-knap target xs)
      (if (null? xs)
          1
          (if (< target (car xs))
    	  (multiply-knap target (cdr xs))
    	  (max (* (car xs)
    		  (multiply-knap (quotient target (car xs)) (cdr xs)))
    	       (multiply-knap target (cdr xs))))))
    (define (nsd-new n)
      (let* ((fs (reverse (factors n)))
    	 (b (multiply-knap (integer-sqrt n) fs))
    	 (a (quotient n b)))
        (list a b)))
    
  4. Paul said

    Here a description for a version of nsd, which is superfast. It solves Euler 266 in Python in 7 seconds! I implemented a C version of an improved “Matthew” algo, which solved it in about 100 minutes (a comparable Python version would have taken 100 hours).

    Split the factors of number n in 2 equal parts, calculate all divisors for the 2 parts and sort both sets of divisors. Loop over both sets of divisors (one ascending and the other descending) and make sure the product is less equal isqrt(n). If the product is too high advance the descending, otherwise the ascending. Keep track of the highest prod below isqrt(n).

  5. Rudolph-san said

    @Paul .. Could you give sample snippet or pseudo code of your idea? Especially about the splitting of divisors?

  6. Paul said

    @Rudolf-san. Here is Python code. You need routines for integer square root (isqrt) and factoring (rho_factors).
    In fact you do not split the divisors, but you split the factors and then calculate all divisors for them. Every divisor for n is a product of two divisors from either part. You are only interested in these products close to isqrt(n).

    def divisors(facs):
        factors = [(1,) + tuple(accumulate(g, mul)) for _, g in groupby(facs)]  # e.g.: [(1, 2), (1, 3), (1, 5) .... (1, 181)]
        div = [1]
        for g in factors:
            div = [d * f for d in div for f in g]
        return div
    
    
    def nsd5(number):
        """ output: nearly square divisor
            method: split factors in 2 equal parts
                    create divisors with small factors and sort (descending)
                    create divisors with large factors (<= ulimit) and sort
                    loop over large divisors and small divisors and search
                       for highest product <= ulimit
        """
        ulimit = isqrt(number)
        facs = rho_factors(number)
        mid = len(facs) // 2
        descending = reversed(sorted(divisors(facs[:mid])))
        ascending = iter(sorted(divisors(facs[mid:])))
        best = 0
        desc = next(descending)
        while True:
            for asc in ascending:
                prod = asc * desc
                if prod > best:
                    if prod <= ulimit:
                        best = prod
                    else:
                        break
            else:
                break  # while
            for desc in descending:
                prod = asc * desc
                if prod <= ulimit:
                    if prod > best:
                        best = prod
                    break
            else:
                break  # while
        return best
    

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: