Phil Harvey’s Puzzle
November 15, 2011
According to the binomial theorem, there are 12870 ways to choose 8 items from a set of 16:
(define (choose n k)
(if (zero? k) 1
(* n (/ k) (choose (- n 1) (- k 1)))))
> (choose 16 8)
12870
A recursive algorithm enumerates the k-subsets of n items. Consider the first element of the list. There are two kinds of subsets: those that include the first element, and those that don’t. Those that do can be enumerated by combining the first element with all (k−1)-subsets of the rest of the elements. Those that don’t can be enumerated by taking all the k-subsets of the rest of the elements. The recursion stops when k is zero or the list is empty. Here is our version of the algorithm:
(define (combs n xs)
(define (comb x xs)
(if (null? xs) xs
(if (pair? (car xs))
(cons (append (list x) (car xs)) (comb x (cdr xs)))
(cons (list x (car xs)) (comb x (cdr xs))))))
(if (or (zero? n) (null? xs)) (list)
(if (= n 1) xs
(append (comb (car xs) (combs (- n 1) (cdr xs)))
(combs n (cdr xs))))))
Now it’s easy. The sum of the numbers from 1 to 16 is 136, the sum of their squares is 1496, and the sum of their cubes is 18496, so we enumerate all 12870 subsets and check for half of each of the three sums:
> (let loop ((xss (combs 8 (range 1 17))))
(let* ((xs (car xss))
(s1 (apply + xs))
(s2 (apply + (map square xs)))
(s3 (apply + (map cube xs))))
(if (and (= s1 68) (= s2 748) (= s3 9248)) xs
(loop (cdr xss)))))
(1 4 6 7 10 11 13 16)
The other partition is thus (2 3 5 8 9 12 14 15).
We used range from the Standard Prelude. The definitions of square and cube are obvious, and are given at http://programmingpraxis.codepad.org/90qQobLz, where you can run the program.
Programming Praxis 
Here’s my Python interactive session, with one (typo, error) removed.
#include <stdio.h> #include <math.h> void process(int s) { int asum, asum2, asum3 ; int bsum, bsum2, bsum3 ; int i ; asum = asum2 = asum3 = 0 ; bsum = bsum2 = bsum3 = 0 ; for (i=0; i<16; i++) { int val, val2, val3 ; val = i ; val2 = val * val ; val3 = val2 * val ; if (s & (1<<i)) { asum += val ; asum2 += val2 ; asum3 += val3 ; } else { bsum += val ; bsum2 += val2 ; bsum3 += val3 ; } } if (asum == bsum && asum2 == bsum2 && asum3 == bsum3) { printf("%d %d %d\n", asum, asum2, asum3) ; printf("A: ") ; for (i=0; i<16; i++) if (s & (1<<i)) printf("%d ", i+1) ; printf("B: ") ; for (i=0; i<16; i++) if ((s & (1<<i)) == 0) printf("%d ", i+1) ; printf("\n") ; } } main() { int s ; for (s=0; s<(1<<15); s++) process(s) ; }If you run my version of the program, you’ll see that there are two unique solutions, only one of which divides { 1..16 } into subsets of equal size. Both solutions have the same sums though. Interesting.
Ah. That teaches me to code before breakfast. I have a small error in the program above: the line which reads val = i should be “val = i + 1 ;”.
When that change is made, it finds the single, unique solution to the problem. My hint that something was wrong was that the sum of the numbers was 60. The sum of the numbers from 1 to 16 should be 136, and since they are divided into two partitions, they should each sum to 68.
Since Michi and kernelbob already posted some Python, here’s a Common Lisp solution:
(defun combinations (xs k) (cond ((or (zerop k) (null xs)) nil) ((= k 1) (mapcar #'list xs)) (t (append (mapcar #'(lambda (c) (cons (first xs) c)) (combinations (rest xs) (1- k))) (combinations (rest xs) k))))) (defun sums (n) (list (/ (* n (1+ n)) 2) (/ (* n (1+ n) (1+ (* 2 n))) 6) (expt (/ (* n (1+ n)) 2) 2))) (defun sum (seq) (reduce #'+ seq :initial-value 0)) (defun square (x) (* x x)) (defun cube (x) (* x x x)) (defun prop (subseq) (equal (mapcar #'(lambda (x) (/ x 2)) (sums 16)) (list (sum subseq) (sum (mapcar #'square subseq)) (sum (mapcar #'cube subseq))))) (let ((xs (loop for i from 1 to 16 collect i))) (dolist (c (remove-if-not #'prop (combinations xs 8))) (print (cons c (list (set-difference xs c))))))and here’s a Haskell one:
import Data.List ((\\)) combinations :: [Integer] -> Integer -> [[Integer]] combinations _ 0 = [] combinations [] _ = [] combinations xs 1 = map (:[]) xs combinations (x:xs) k = map (x:) (combinations xs (k - 1)) ++ combinations xs k sums :: Integer -> [Integer] sums n = [(n * (n + 1)) `div` 2, (n * (n + 1) * (2 * n + 1)) `div` 6, (n * (n + 1) `div` 2) ^2] prop :: [Integer] -> Bool prop xs = map (`div` 2) (sums 16) == map sum [xs, map (^2) xs, map (^3) xs] main :: IO () main = mapM_ printPairs filteredCombs where printPairs = print . (\ c -> (c, [1..16] \\ c)) filteredCombs = filter prop $ combinations [1..16] 8In both, I calculated the sums of the integers, squares of integers, and cubes of integers in {1, …, 16} via formulae for the first three types of figurate numbers.
Here’s a ruby version …
def sum(x) x.inject(0) { |sum, v| sum += v } end def sum_squares(x) x.inject(0) { |sum, v| sum += v**2 } end def sum_cubes(x) x.inject(0) { |sum, v| sum += v**3 } end x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] sum_x = sum(x) sum_squares_x = sum_squares(x) sum_cubes_x = sum_cubes(x) x.permutation(8) do |p| if (sum(p) == sum_x/2) && (sum_squares(p) == sum_squares_x/2) && (sum_cubes(p) == sum_cubes_x/2) puts "answer = #{p.sort} #{x-p}" break; end endThis problem can be formulated as a (binary) integer programming problem.
For example, using GMPL
var x{i in 1..16}, binary; s.t. pow{j in 0..3}: sum{i in 1..16} i^j * x[i] = (sum{i in 1..16} i^j)/2; solve; display{i in 1..16} x[i]; end;I just found the “snoob” function over here. http://www.steike.com/code/useless/evil-c/ That led to this somewhat golfed C solution.
Basically, snoob(n) returns the next integer after n with the same number of bits set. If we use a bit mask to denote the members of the first partition, we can use snoob() to iterate through the partitions directly.
for (i = 0xFF; i <= 0xFF00; i = snoob(i)) { …}
// 0xFF is the first partition, {1, 2, 3, 4, 5, 6, 7, 8}.
// 0xFF00 is the last partition, {9, 10, 11, 12, 13, 14, 15, 16}.
snoob() should be blazingly fast since it's about 7 arithmetic operations and no memory references.
#include <stdio.h>
#include <string.h>
unsigned snoob(unsigned x) {
unsigned smallest = x & -x;
unsigned ripple = x + smallest;
unsigned ones = x ^ ripple;
ones = (ones >> 2) / smallest;
return ripple | ones;
}
int main()
{
unsigned i, j, k, l;
for (i = 0xFF; i <= 0xFF00; i = snoob(i)) {
unsigned sums[2][3] = { { 0, 0, 0 }, { 0, 0, 0 } };
for (j = 0; j < 16; j++)
for (k = 0, l = 1; k < 3; k++)
sums[!(i & 1 << j)][k] += l *= j + 1;
if (!memcmp(sums[0], sums[1], sizeof sums[0])) {
for (j = 0; j < 16; j++)
if (i & 1 << j)
printf("%d ", j + 1);
printf("\n");
}
}
return 0;
}
Oops. Formatted.
#include <stdio.h> #include <string.h> unsigned snoob(unsigned x) { unsigned smallest = x & -x; unsigned ripple = x + smallest; unsigned ones = x ^ ripple; ones = (ones >> 2) / smallest; return ripple | ones; } int main() { unsigned i, j, k, l; for (i = 0xFF; i <= 0xFF00; i = snoob(i)) { unsigned sums[2][3] = { { 0, 0, 0 }, { 0, 0, 0 } }; for (j = 0; j < 16; j++) for (k = 0, l = 1; k < 3; k++) sums[!(i & 1 << j)][k] += l *= j + 1; if (!memcmp(sums[0], sums[1], sizeof sums[0])) { for (j = 0; j < 16; j++) if (i & 1 << j) printf("%d ", j + 1); printf("\n"); } } return 0; }A haskell one-liner in ghci:
And a little sorter solution with list comprehension:
@neez Nice Haskell solution. You can also do it without precalculating the constants.
@Jan Van lent
Wow, beautiful soution, especially using the i=0 case to partition into equal sized sets!
@neez Well, it is only a small change from your solution. I used the same trick in the GMPL solution. Note that the constraint that the sets be of equal size is actually redundant. You can try with [1..3], [2..3] and even just 3 and they all give the same solution.
Create 2 (random) partitions, and swap between them, making sure the difference is less. Keep doing this until we get a difference of 0.
import random gen = random.Random() def diff(s, t): return sum(abs(sum((k**x for k in s)) - sum((i**x for i in t))) for x \ in range(1,4)) def swap(s, t): start = diff(s,t) end = start + 1 while end > start: ind1 = gen.randint(0, 7) ind2 = gen.randint(0, 7) s[ind1], t[ind2] = t[ind2], s[ind1] end = diff(s,t) def main(): x = [i for i in range(1,17)] s = random.sample(x, 8) t = [k for k in x if k not in s] while diff(s,t) != 0: swap(s,t) s.sort() t.sort() print(s) print(t) if __name__ == '__main__': main()public class HarveyPuzzle { public static final int ASIZE=16; public HarveyPuzzle(){ puzzleStart(); } public static void main(String[] args) { new HarveyPuzzle(); } private void puzzleStart() { boolean [] pick = new boolean[ASIZE]; recursiveFinder(pick, 0); } private void recursiveFinder(boolean [] pick, int k){ if(k==pick.length-1){ isSolution(pick); pick[k]=!pick[k]; isSolution(pick); } else { recursiveFinder(pick, k+1); pick[k]=!pick[k]; recursiveFinder(pick, k+1); } } private void isSolution(boolean [] pick){ if(sum(pick, true)==sum(pick, false) && squaresSum(pick, true)==squaresSum(pick, false) && cubesSum(pick, true)==cubesSum(pick, false)){ System.out.println("Solution found!"); printSolution(pick); } } private int sum(boolean [] pick, boolean w) { int somma=0; for(int i=0;i<pick.length;i++){ if(pick[i]==w) somma+=i+1; } return somma; } private int squaresSum(boolean [] pick, boolean w) { int somma=0; for(int i=0;i<pick.length;i++){ if(pick[i]==w) somma+=Math.pow(i+1, 2); } return somma; } private int cubesSum(boolean [] pick, boolean w) { int somma=0; for(int i=0;i<pick.length;i++){ if(pick[i]==w) somma+=Math.pow(i+1, 3); } return somma; } private void printSolution(boolean[] pick) { System.out.print("["); for(int i=0;i<pick.length;i++){ if(pick[i]){ System.out.print((i+1)+" "); } } System.out.print("], ["); for(int i=0;i<pick.length;i++){ if(!pick[i]){ System.out.print((i+1)+" "); } } System.out.println("]\n"); System.err.println(sum(pick, true) + "==" + sum(pick, false)); System.err.println(squaresSum(pick, true) + "==" + squaresSum(pick, false)); System.err.println(cubesSum(pick, true) + "==" + cubesSum(pick, false)+"\n\n"); } }package com.algorithm;
import java.util.Arrays;
public class SameSumSquareCube {
private static int[] set1;
private static int[] set2;
public static void main(String[] args) {
set1 = new int[] {1, 2, 3, 4, 5, 6, 7, 8};
set2 = new int[] {9, 10, 11, 12, 13, 14, 15, 16};
recursive(7);
System.out.println(“completed”);
}
private static void recursive(int count) {
for (int i=count; i>=0; i–) {
for (int j=7; j>=0; j–) {
int temp = set1[i];
set1[i] = set2[j];
set2[j] = temp;
if (isCubeSumEqual(set1, set2) && isSquareSumEqual(set1, set2) && isSumEqual(set1, set2)) {
System.out.println(Arrays.toString(set1));
System.out.println(Arrays.toString(set2));
}
}
recursive(–count);
}
}
private static boolean isCubeSumEqual(int[] set1, int[] set2) {
int set1CubeSum = 0;
for (int number : set1) {
set1CubeSum += number * number * number;
}
int set2CubeSum = 0;
for (int number : set2) {
set2CubeSum += number * number * number;
}
return set1CubeSum == set2CubeSum;
}
private static boolean isSquareSumEqual(int[] set1, int[] set2) {
int set1SquareSum = 0;
for (int number : set1) {
set1SquareSum += number * number;
}
int set2SquareSum = 0;
for (int number : set2) {
set2SquareSum += number * number;
}
return set1SquareSum == set2SquareSum;
}
private static boolean isSumEqual(int[] set1, int[] set2) {
int set1Sum = 0;
for (int number : set1) {
set1Sum += number;
}
int set2Sum = 0;
for (int number : set2) {
set2Sum += number;
}
return set1Sum == set2Sum;
}
}
[…] the weekend, I was amusing myself with this problem from the always entertaining Programming Praxis. The problem is to partition the set {1, 2, […]