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.

About these ads

Pages: 1 2

18 Responses to “Phil Harvey’s Puzzle”

  1. kernelbob said

    Here’s my Python interactive session, with one (typo, error) removed.

    ~> python
    Python 2.7.1+ (r271:86832, Apr 11 2011, 18:13:53) 
    [GCC 4.5.2] on linux2
    Type "help", "copyright", "credits" or "license" for more information.
    >>> a = range(1, 17)
    >>> a
    [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    >>> from itertools import combinations
    >>> s1 = sum
    >>> s2 = lambda s: sum(i**2 for i in s)
    >>> s3 = lambda s: sum(i**3 for i in s)
    >>> s3(a)
    18496
    >>> s2(a)
    1496
    >>> s1(a)
    136
    >>> for p in combinations(a, 8):
    ...  if all(f(p) == t / 2 for f, t in zip((s1, s2, s3), (136, 1496, 18496))):
    ...   print p
    ... 
    (1, 4, 6, 7, 10, 11, 13, 16)
    (2, 3, 5, 8, 9, 12, 14, 15)
    >>> 
    
  2. 
    #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) ;
    }
    
  3. 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.

  4. 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.

  5. Graham said

    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] 8
    

    In 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.

  6. slabounty said

    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 
    end
    
  7. Jan Van lent said

    This 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;
    
  8. kernelbob said

    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;
    }

  9. kernelbob said

    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;
    }
    
  10. neez said

    A haskell one-liner in ghci:

    Prelude Data.List> let ps p = sum.map (^p) in filter (\x -> and [8==length x,68==ps 1 x,748==ps 2 x,9248==ps 3 x]) $ subsequences [1..16]
    [[2,3,5,8,9,12,14,15],[1,4,6,7,10,11,13,16]]
    
  11. neez said

    And a little sorter solution with list comprehension:

    Prelude Data.List> [x | x <- subsequences [1..16], let ps p = sum.map (^p), 8==length x, 68==ps 1 x, 748==ps 2 x, 9248==ps 3 x]
    [[2,3,5,8,9,12,14,15],[1,4,6,7,10,11,13,16]]
    
  12. Jan Van lent said

    @neez Nice Haskell solution. You can also do it without precalculating the constants.

    Prelude Data.List> let ps p = sum.map (^p) in filter (\x -> and [ps i [1..16]==2*(ps i x) | i <- [0..3]]) $ subsequences [1..16]
    [[2,3,5,8,9,12,14,15],[1,4,6,7,10,11,13,16]]
    
  13. neez said

    @Jan Van lent
    Wow, beautiful soution, especially using the i=0 case to partition into equal sized sets!

  14. Jan Van lent said

    @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.

  15. Yuushi said

    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()
    
  16. Simone said
    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");
    	}
    	
    }
    
  17. Vaibhav Bansal said

    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;
    }
    }

  18. […] the weekend, I was amusing myself with this problem from the always entertaining Programming Praxis. The problem is to partition the set {1, 2, […]

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

Follow

Get every new post delivered to your Inbox.

Join 643 other followers

%d bloggers like this: