Three-Way Minimum Sum Partitions
November 22, 2019
Our function f takes a single argument n and returns a list of three elements (x y z) in increasing order; for instance, (f 6) returns (1 2 3). If n is prime, the function returns (1 1 n); if n = p * q is semi-prime, the function returns (1 p q); and if n has three prime factors, the function returns (x y z). Life gets harder when n has four or more factors.
Our first function uses brute force to compute all possible combinations of x, y and z and find the minimal sum:
(define (f1 n) (let ((xyz (list n n n))) (do ((xs (divisors n) (cdr xs))) ((null? xs) (sort < xyz)) (do ((ys (divisors (/ n (car xs))) (cdr ys))) ((null? ys)) (let ((z (/ n (car xs) (car ys)))) (when (< (+ (car xs) (car ys) z) (sum xyz)) (set! xyz (list (car xs) (car ys) z))))))))
That works fine for n with 1, 2, or 3 factors, or more, though it takes a while when n has lots of prime factors (that last number has 13,440 divisors):
> (f1 17) (1 1 17) > (f1 13290059) (1 3119 4261) > (f1 (* 5 17 83)) (5 17 83) > (f1 (* 2 3 5 83)) (5 6 83) > (time (f1 (* 2 2 2 2 2 2 3 3 3 3 5 5 5 7 7 11 13 17 19 23))) (time (f1 (* 2 ...))) 49 collections 7.661519834s elapsed cpu time, including 0.003838511s collecting 7.686354068s elapsed real time, including 0.004012308s collecting 412089984 bytes allocated, including 412857248 bytes reclaimed (32292 32300 32340)
Our second function finds the largest divisor less than or equal to the cube root of n, then the largest divisor less than or equal to the square root of the remaining cofactor of n, and then whatever is left:
(define (f2 n)
(let* ((x (maxle (iroot 3 n) (divisors n)))
(y (maxle (iroot 2 (/ n x)) (divisors (/ n x))))
(z (/ n x y)))
(sort < (list x y z))))
> (f2 17) (1 1 17) > (f2 13290059) (1 3119 4261) > (f2 (* 5 17 83)) (5 17 83) > (f2 (* 2 3 5 83)) (3 10 83) > (time (f2 (* 2 2 2 2 2 2 3 3 3 3 5 5 5 7 7 11 13 17 19 23))) (time (f2 (* 2 ...))) no collections 0.110332680s elapsed cpu time 0.110471207s elapsed real time 1127552 bytes allocated (32292 32300 32340)
That’s faster, but it doesn’t always work, as the example of 2 * 3 * 5 * 83 = 2490 shows; the large factor overwhelms the rest of the calculation. To try to overcome large factors, my third function makes two calculations — the f2 calculation, plus a new calculation where I first remove the largest factor then take the two divisors on either side of the square root of the remaining cofactor — and compares the two, taking the smaller:
(define (f3 n) (let* ((x1 (maxle (iroot 3 n) (divisors n))) (y1 (maxle (iroot 2 (/ n x1)) (divisors (/ n x1)))) (z1 (/ n x1 y1)) (x2 (apply max (factors n))) (y2 (maxle (iroot 2 (/ n x2)) (divisors (/ n x2)))) (z2 (/ n x2 y2))) (if (< (+ x1 y1 z1) (+ x2 y2 z2)) (sort < (list x1 y1 z1)) (sort < (list x2 y2 z2)))))
But that still doesn’t work. Consider n = 1890 = 2 * 3 * 3 * 3 * 5 * 7: (f1 1890) => (9 14 15), a sum of 38; (f2 1890) => (9 10 21), a sum of 40; and (f3 1890) => (7 15 18), a different sum of 40.
So I’m not sure what to do next. Relying on the brute force calculation is unsatisfying, and can be very slow, but nothing better comes to mind. You can run the program at https://ideone.com/wV0vHJ.
I wonder if choosing a series of random partitions might work….
Programming Praxis 
The method below is much faster than brute force. First determine all divisors and sort them in decreaing order. Loop for x starting from the cubic root downwards. Divide n by x, determine the sqrt of the result and check if it possible to improve the sum. Break if not possible. Otherwise loop downwards from this sqrt and find the first possible solution.
I will post full code later on ideone.
def threeway3(n): divs = sorted(divisors(n), reverse=True) n3 = iroot(3, n) min_val, min_sol = inf, None for x in divs: if x > n3: continue m = n // x sqrtm = isqrt(m) if x + 2 * sqrtm > min_val: break for y in divs: if y < x: break if y > sqrtm: continue if m % y == 0: z = m // y if x + y + z < min_val: min_val = x + y + z min_sol = (x, y, z) break return min_sol, min_valFull code with some changes is on Ideone.
@Paul: nice – I was thinking on similar lines, but didn’t think of the cutoff on line 10. I think that break on line 22 can be outdented (any smaller divisor than y will have a larger sum). It would be nice to do a binary search through the divisors to find the start points of the iterations.
@matthew. I saw also that I could outdent that line (see my Ideone post)
@Paul: So I see, good stuff. I’m wondering if we can remove cubic factors from n before finding the minimum sum, ie. if a+b+c = n is minimal, then is it the case that ka+kb+kc = k^3n is minimal?
@matthew. An interesting idea, but unfortunately there is a counterexample. Take 24=83. The best solution is x,y,z = 2, 3, 4 with sum 9. If you solve for n=3 with solution 1,1,3 (sum 5) the total sum is 25=10.
@matthew. Formatting caught me here, Of course, I meant 24=8×3 and 2×5=10.
@Paul: nice counterexample. Probably wouldn’t have helped much in the general case anyway.
Here’s a recursive solution in Python that supports an arbitrary number of factors, as opposed to 3 specified in the problem.
def prime_factorize(n): assert n >= 2 factors = [] x = 2 while x * x <= n: if n % x == 0: factors.append(x) n //= x continue x += 1 + (x & 1) factors.append(n) return tuple(factors) def calc_divisors(n): divisors = {1} for factor in prime_factorize(n): divisors.update([factor * d for d in divisors]) return tuple(divisors) def min_sum_factorize(n, num_factors): if num_factors == 1: return n, best = (n,) + (1,) * (num_factors - 1) best_sum = sum(best) for x in calc_divisors(n): if pow(x, num_factors) > n: continue if n % x == 0: rest = min_sum_factorize(n // x, num_factors - 1) total = sum(rest) + x if total < best_sum: best = (x,) + rest best_sum = total return best print(min_sum_factorize(2490, 3)) print(min_sum_factorize(1890, 3)) print(min_sum_factorize(33731641944000, 3)) print(min_sum_factorize(2490, 4)) print(min_sum_factorize(1890, 4)) print(min_sum_factorize(33731641944000, 4))Output:
Here’s the recursive solution I posted above, with @Paul’s short-circuiting optimization. It runs noticeably faster now :-), even with the naive integer root function I wrote.
def iroot(n, k): """A naive algorithm for calculating the degree k integer root of n""" x = 1 while pow(x, k) < n: x *= 2 while pow(x, k) > n: x -= 1 return x def prime_factorize(n): assert n >= 2 factors = [] x = 2 while x * x <= n: if n % x == 0: factors.append(x) n //= x continue x += 1 + (x & 1) factors.append(n) return tuple(factors) def calc_divisors(n): divisors = {1} for factor in prime_factorize(n): divisors.update([factor * d for d in divisors]) return tuple(divisors) def min_sum_factorize(n, num_factors): num_factors_minus_one = num_factors - 1 if num_factors == 1: return n, best = (n,) + (1,) * num_factors_minus_one best_sum = sum(best) for x in sorted(calc_divisors(n), reverse=True): if pow(x, num_factors) > n: continue if x + num_factors_minus_one * iroot(n // x, num_factors_minus_one) > best_sum: break if n % x == 0: rest = min_sum_factorize(n // x, num_factors_minus_one) total = sum(rest) + x if total < best_sum: best = (x,) + rest best_sum = total return best for n in (2490, 1890, 33731641944000): for num_factors in (3, 4, 5): result = min_sum_factorize(n, num_factors) print(f'min_sum_factorize({n}, {num_factors}) = {result}')Output:
Klong solution
In addition to the brute force approach, I tried the following 2 alternatives:
1) Iterated the first loop from 1 to the floored cube root of the initial number; and
2) Iterated the first loop from 1 to the floored cube root of the initial number and iterated the 2nd loop from 1 to the floored square root of the quotient of the initial number and the 1st iterator.
The run shows one of the combinations which yields one of the least sums as well as the number of combinations before finding the minimum.
go::{[l]; l::l@<{+/x}'l::{x}'go1(x); (#l),,*l} go1::{[l n q]; l::[]; n::x;n{:[q=_q::n%x; l::go2(l; x; _q);""]; x+1}:*1; l} go2::{[l2 y2 n2 q2]; l2::x; y2::y; n2::z; n2{:[q2=_q2::n2%x; l2::l2,,y2,x,_q2; 1]; x+1}:*1; l2} go3::{[l]; l::l@<{+/x}'l::{x}'go4(x); (#l),,*l} go4::{[l n q]; l::[]; n::x; (_n^1%3){:[q=_q::n%x; l::go5(l; x; _q);""]; x+1}:*1; l} go5::{[l2 y2 n2 q2]; l2::x; y2::y; n2::z; n2{:[q2=_q2::n2%x; l2::l2,,y2,x,_q2; 1];x+1}:*1; l2} go6::{[l]; l::l@<{+/x}'l::{x}'go7(x); (#l),,*l} go7::{[l n q]; l::[]; n::x; (_n^1%3){:[q=_q::n%x; l::go8(l; x; _q);""]; x+1}:*1; l} go8::{[l2 y2 n2 q2]; l2::x; y2::y; n2::z; (_n2^1%2){:[q2=_q2::n2%x; l2::l2,,y2,x,_q2; 1]; x+1}:*1; l2} n::#l::go(1890); .p("Brute force: Answer - [",(l@1),"], # Entries - ",,(l@0)); .p("") n::#l::go3(1890); .p("Square root: Answer - [",(l@1),"], # Entries - ",,l@0); .p("") n::#l::go6(1890); .p("Cube/Square roots: Answer - [",(l@1),"], # Entries - ",,l@0); .p("")Results from the Klong run:
[Brute force: Answer – [ 9 14 15 ], # Entries – 270]
[Square root: Answer – [ 9 15 14 ], # Entries – 140]
[Cube/Square roots: Answer – [ 9 14 15 ], # Entries – 70]