Programs that Write Programs

Ever heard of the Collatz Conjecture? Neither had I until about two hours (of study time I will never get back) ago.

Wasting time yesterday, I came across it in an entry on Roshan's blog. Of course, trolling the internets is a dumb thing to do when you have a crushing final approaching - and this was a very good lesson in why I SHOULDN'T take breaks in this way. I intended to pick this up again after the test, but of course I couldn't get it out of my head.

See, there's this function that seems to return 1 no matter what (positive?) integer you feed it.

It works like this: if the number you feed it is 1, return 1 - you're done! If the number you feed it is even, return the result of applying the function to half that number. And if the number is odd, multiply it by 3, add 1, and then apply the function.

Here's a C-like program (stolen from Roshan's blog) that implements it:


int collatz(int n)
{
if(n == 1)
  return n;
else if (n % 2 == 0)
  return collatz(n/2);
else
  return collatz(n*3+1);
}


Now, Roshan and Will Byrd spent a long time talking about this one day and came up with (what looks to be) another such function.

And here's Roshan's version in Scheme:

(define C
  (lambda (n)
    (cond
      ((= n 1) 1)
      ((= 0 (modulo n 3)) (C (quotient n 3)))
      ((= 1 (modulo n 3)) (C (- (* n 4) 1)))
      (else (C (+ (* n 4) 1))))))

If you can read the C code above, this Scheme code should be pretty easy to follow too. The idea here is that the "divisor" number is 3 this time instead of 2, and there are now two conditions rather than one. So, if the number is 1, it returns 1 and you're finished. If the number is evenly divisible by three, then return the result of the function applied to the number divided by 3. If the number divides by 3 with a remainder of 1, then return the function applied to one less than the number times 4. ELSE: return the result of the function applied to one more than the number times 4.

Roshan then found YET ANOTHER such function later on:

(define C5
  (lambda (n)
    (cond
      [(= n 1) 1]
      [(= 0 (modulo n 5)) (C5 (quotient n 5))]
      [(= 1 (modulo n 5)) (C5 (+ (* 6 n) 4))]
      [(= 2 (modulo n 5)) (C5 (+ (* 6 n) 3))]
      [(= 3 (modulo n 5)) (C5 (+ (* 6 n) 2))]
      [else (C5 (+ (* 6 n) 1))])))

(NOTE: this is actually my Scheme code, but the idea was taken from Roshan's blog entry linked above.)

At the end of his entry, Roshan speculates that there might be an infinite number of such functions. The interesting thing about them is that no one has been able to verify that this property (of returning 1 for whatever integer you feed it) actually holds for all integers. That's why it's still a "conjecture." It seems intuitively right - that is, if you work out a number of examples, you get the feeling it's probably true - but feelings aren't the same thing as a real proof. In fact, there are several historical examples of things like this that seemed to be true until people started testing really huge numbers - and then they strangely weren't true anymore.

At any rate, no one has come up with a proof for the original function yet. Nor, it goes without saying, does anyone have a proof for Roshan's two functions. But Roshan tested them for all integers up to 100000 on the computer, and they do indeed work for that range.

Alright, so I wanted to play too - and I took a good hard look at Roshan's two functions trying to come up with a generalization. I did come up with one that seemed to work - but unfortunately it broke down for some of the examples.

But I'm pleased to say that I DID discover one of these functions of my own:

(define C7
  (lambda (n)
   (cond
    [(= n 1) 1]
    [(= 0 (modulo n 7)) (C7 (quotient n 7))]
    [(= 1 (modulo n 7)) (C7 (+ (* 8 n) 6))]
    [(= 2 (modulo n 7)) (C7 (+ (* 8 n) 5))]
    [(= 3 (modulo n 7)) (C7 (+ (* 8 n) 4))]
    [(= 4 (modulo n 7)) (C7 (+ (* 8 n) 3))]
    [(= 5 (modulo n 7)) (C7 (+ (* 8 n) 2))]
    [else (C7 (+ (* 8 n) 1))])))

Now you see why I included my own Scheme code for the "5" example! :-) It should be easy to see what "generalization" I thought I'd come up with. Or rather, the generalization that I thought I'd come up with for odd numbers. It seemed that what was going on was that there was a "divisor" number that characterized the function. (2 in the original example, 3 in Roshan and Will's example and 5 in Roshan's second example.) If the number divides evenly by that number, you recur on the result of dividing the input number by the characterisitic divisor. Otherwise, if the remainder of dividing the input number by the characteristic divisor is at least 2 smaller than the characteristic divisor, then you recur on whatever the difference between the remainder and the divisor is plus the input number multiplied by one more than the characteristic divisor. (<- It's really not possible to make this sentence any clearer; you'll probably have more luck just looking at the code. You'll get the hang of it if you add the first number in each line that starts with a number to the last number in the same line. Notice that they always total whatever the number in the (modulo n x) statement is.)

Alright, well, I thought I'd found a pattern on the odd numbers, so I wrote the following generator program to "make" these functions:

(define C-gen
  (lambda (c n)
   (cond
    [(= n 1) 1]
    [(= 0 (modulo n c)) (C-gen c (quotient n c))]
    [(> (- c (modulo n c)) 1) (C-gen c (+ (* (+ c 1) n) (- c (modulo n c))))]
    [else (C-gen c (+ (* (+ c 1) n) 1))])))

And this does indeed "make" Collatz functions.

Of course, you can make it more "functional" (i.e. so that it returns a function - that is, is a program that writes a program) by "stacking" the lambda arguments. By way of example, here's one that I wrote to do Roshan's test. First you feed it a bound (100000, to keep up with Roshan), and it returns you a function that tests, for a given divisor, whether the function generated according to my generalization for that divisor really does return 1 for all numbers from 1 to 100000:

(define verify-collatz
  (lambda (c bound)
   (cond
    [(zero? bound) #t]
    [(= 1 (C-gen c bound)) (verify-collatz c (- bound 1))]
    [else #f])))

Here's where I had my first heartbreak. It turns out that it doesn't work for 9!

Here are some counterexamples (actually, I just got impatient after 4 minutes or so - maybe they just take a really long time to return): 99991, 99957, 99860, 99841, 99310, 99291, 99277. These don't seem to return 1 for characteristic divisor 9.

Here's a program that prints out "Works for c: 9, n: 99999" for each number that the function returns 1 for:

(define vc-gen
  (lambda (bound)
   (lambda (c)
    (cond
    [(zero? bound) #t]
    [(= 1 (C-gen c bound)) (begin (printf "Works for c: ~s, n: ~s ~%" c bound) ((vc-gen (- bound 1)) c))]
    [else #f]))))

This helps to spot numbers that it "suspiciously gets hung up on." Again, there's no rigorous proof for these numbers, but what it does for 9 on my machine is that it shows that it works for 100000, 99999, 99998 etc. almost immediately, then it just never seems to return for 99991. So I'm making the assumption that it diverges for all the numbers listed above - but again, I don't really know that.

Well, so much for my odd idea. It was a dumb idea anyway. After all - the original function is for 2, and 2 is not an odd number!!!

Give it up - it doesn't work for the primes either. 17 is the counterexample. It cuts out on you at 99992.

However - I did come up with an idea that I haven't found a counterexample to yet: the Fibonacci numbers!!!

I've run tests for the whole series up to 987. 987 is running at the time of this writing. It hasn't returned yet, but I have it in display mode, and it's already in the 41000s without getting hung up. I will do a more rigorous test on a large portion of the series later - but this seems to be the answer.

So, as a VERY TENTATIVE answer to Roshan's question: the answer seems to be "Yes, there is an infinite number of such functions. You can generate such a function for each of the numbers in the fibonacci sequence."

Not that any of this is helping me pass my Algorithms class...

[UPDATE: It just returned true for 987!!! I think I'm on to something.]