Nearly Square Divisors, Revisited

August 23, 2016

The knap function takes a target, which is the logarithm of the square root of n, and a list xs of the logarithms of the factors of n, with multiplicity, in descending order, and returns the logarithm of the largest divisor of n less than the logarithm of the square root. The function is recursive, with the target decreasing and logarithms of the factors being removed from the list. The base case occurs when the list of logarithms of the factors of n is empty, when the function returns 0. If the largest logarithm is larger than the target, the recursion continues with the remaining logarithms. Otherwise, the algorithm returns the maximum of two values: the sum of the largest logarithm and the reduced knapsack problem with the largest logarithm subtracted from the target and removed from the list, or the recursion with the remaining logarithms. The code is simpler than the explanation:

(define (knap target xs)
  (if (null? xs) 0
    (if (< target (car xs))
        (knap target (cdr xs))
        (max (+ (car xs) (knap (- target (car xs)) (cdr xs)))
             (knap target (cdr xs))))))

Then the nearly square divisors are computed by factoring n, arranging the logarithms in descending order, calculating the target as half the sum of the logarithms (which is the logarithm of the square root of n), calling the knap function, taking the anti-logarithm, and finally dividing to find the other divisor:

(define (nsd n)
  (let* ((fs (map log (reverse (factors n))))
         (b (inexact->exact (round (exp
              (knap (/ (sum fs) 2) fs)))))
         (a (/ n b)))
    (list a b)))

Here’s the example from the previous exercise:

> (time (nsd 224403121196654400))
(time (nsd 224403121196654400))
    8 collections
    0.121862989s elapsed cpu time, including 0.000410511s collecting
    0.124334069s elapsed real time, including 0.000441672s collecting
    65259936 bytes allocated, including 67360576 bytes reclaimed
(473753280 473670855)

This sadly runs into a problem; the logarithms lose precision, and if n is too large, the program fails to find the correct solution. It will take something like the unlimited-precision floating-point functions in gmp to make it work for very large numbers.

You can run the program at http://ideone.com/n5Wcja, where you will also see the sum and factors functions.

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 )

Google photo

You are commenting using your Google 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 )

Connecting to %s

%d bloggers like this: