Hailstones

February 17, 2012

Here is our simple function producing the hailstone sequence:

(define (hailstone n)
  (let loop ((n n) (c '()))
    (if (= n 1)
        (reverse (cons 1 c))
        (loop (if (even? n)
                  (quotient n 2)
                  (+ (* 3 n) 1))
              (cons n c)))))

For example:

> (hailstone 13)
(13 40 20 10 5 16 8 4 2 1)

The longest hailstone sequence for x0 < 1000000 is for a chain starting at 837799 which has 525 elements and reaches a maximum of 2974984576 at its 59th element (marked in bold italic on the tenth line):

> (hailstone 837799)
(837799 2513398 1256699 3770098 1885049 5655148 2827574
1413787 4241362 2120681 6362044 3181022 1590511 4771534
2385767 7157302 3578651 10735954 5367977 16103932 8051966
4025983 12077950 6038975 18116926 9058463 27175390 13587695
40763086 20381543 61144630 30572315 91716946 45858473
137575420 68787710 34393855 103181566 51590783 154772350
77386175 232158526 116079263 348237790 174118895 522356686
261178343 783535030 391767515 1175302546 587651273
1762953820 881476910 440738455 1322215366 661107683
1983323050 991661525 2974984576 1487492288 743746144
371873072 185936536 92968268 46484134 23242067 69726202
34863101 104589304 52294652 26147326 13073663 39220990
19610495 58831486 29415743 88247230 44123615 132370846
66185423 198556270 99278135 297834406 148917203 446751610
223375805 670127416 335063708 167531854 83765927 251297782
125648891 376946674 188473337 565420012 282710006 141355003
424065010 212032505 636097516 318048758 159024379 477073138
238536569 715609708 357804854 178902427 536707282 268353641
805060924 402530462 201265231 603795694 301897847 905693542
452846771 1358540314 679270157 2037810472 1018905236
509452618 254726309 764178928 382089464 191044732 95522366
47761183 143283550 71641775 214925326 107462663 322387990
161193995 483581986 241790993 725372980 362686490 181343245
544029736 272014868 136007434 68003717 204011152 102005576
51002788 25501394 12750697 38252092 19126046 9563023
28689070 14344535 43033606 21516803 64550410 32275205
96825616 48412808 24206404 12103202 6051601 18154804 9077402
4538701 13616104 6808052 3404026 1702013 5106040 2553020
1276510 638255 1914766 957383 2872150 1436075 4308226
2154113 6462340 3231170 1615585 4846756 2423378 1211689
3635068 1817534 908767 2726302 1363151 4089454 2044727
6134182 3067091 9201274 4600637 13801912 6900956 3450478
1725239 5175718 2587859 7763578 3881789 11645368 5822684
2911342 1455671 4367014 2183507 6550522 3275261 9825784
4912892 2456446 1228223 3684670 1842335 5527006 2763503
8290510 4145255 12435766 6217883 18653650 9326825 27980476
13990238 6995119 20985358 10492679 31478038 15739019
47217058 23608529 70825588 35412794 17706397 53119192
26559596 13279798 6639899 19919698 9959849 29879548 14939774
7469887 22409662 11204831 33614494 16807247 50421742
25210871 75632614 37816307 113448922 56724461 170173384
85086692 42543346 21271673 63815020 31907510 15953755
47861266 23930633 71791900 35895950 17947975 53843926
26921963 80765890 40382945 121148836 60574418 30287209
90861628 45430814 22715407 68146222 34073111 102219334
51109667 153329002 76664501 229993504 114996752 57498376
28749188 14374594 7187297 21561892 10780946 5390473 16171420
8085710 4042855 12128566 6064283 18192850 9096425 27289276
13644638 6822319 20466958 10233479 30700438 15350219
46050658 23025329 69075988 34537994 17268997 51806992
25903496 12951748 6475874 3237937 9713812 4856906 2428453
7285360 3642680 1821340 910670 455335 1366006 683003 2049010
1024505 3073516 1536758 768379 2305138 1152569 3457708
1728854 864427 2593282 1296641 3889924 1944962 972481
2917444 1458722 729361 2188084 1094042 547021 1641064 820532
410266 205133 615400 307700 153850 76925 230776 115388 57694
28847 86542 43271 129814 64907 194722 97361 292084 146042
73021 219064 109532 54766 27383 82150 41075 123226 61613
184840 92420 46210 23105 69316 34658 17329 51988 25994 12997
38992 19496 9748 4874 2437 7312 3656 1828 914 457 1372 686
343 1030 515 1546 773 2320 1160 580 290 145 436 218 109 328
164 82 41 124 62 31 94 47 142 71 214 107 322 161 484 242 121
364 182 91 274 137 412 206 103 310 155 466 233 700 350 175
526 263 790 395 1186 593 1780 890 445 1336 668 334 167 502
251 754 377 1132 566 283 850 425 1276 638 319 958 479 1438
719 2158 1079 3238 1619 4858 2429 7288 3644 1822 911 2734
1367 4102 2051 6154 3077 9232 4616 2308 1154 577 1732 866
433 1300 650 325 976 488 244 122 61 184 92 46 23 70 35 106
53 160 80 40 20 10 5 16 8 4 2 1)

We could have calculated the length of each sequence in a loop, but it is better to store the sequence lengths in a cache, build the cache in order, and check at each step if the current sequence is in the cache; that saves an order of magnitude in run time:

(define (longest-sequence n)
  (let ((cs (make-vector 1000000 0))
        (max-x 1) (max-len 1))
    (vector-set! cs 1 1)
    (do ((k 2 (+ k 1)))
        ((= k n) max-x)
        (vector-set! cs k
          (do ((c 0 (+ c 1))
               (s k (if (even? s)
                        (/ s 2)
                        (+ (* 3 s) 1))))
              ((< s k)
                (+ (vector-ref cs s) c))))
        (if (< max-len (vector-ref cs k))
            (begin (set! max-x k)
                   (set! max-len
                         (vector-ref cs k)))))))

And here is the time comparison:

> (time (longest-sequence 1000000))
(time (longest-sequence 1000000))
    21 collections
    1156 ms elapsed cpu time, including 0 ms collecting
    2000 ms elapsed real time, including 0 ms collecting
    92089824 bytes allocated, including 86367800 bytes reclaimed
837799
> (time (apply max (list-of (length (hailstone n)) (n range 1 1000000))))
(time (apply max ...))
    521 collections
    16437 ms elapsed cpu time, including 157 ms collecting
    16000 ms elapsed real time, including 0 ms collecting
    2207238296 bytes allocated, including 2200258672 bytes reclaimed
525

There’s much more interesting recreational mathematics in the hailstone sequence than we have shown here; see Wikipedia, MathWorld or the OEIS for more. You can run our programs at http://programmingpraxis.codepad.org/MUTonE0Q.

About these ads

Pages: 1 2

20 Responses to “Hailstones”

  1. DGel said

    Trying to learn haskell, this seemed like it fits the language very well =) Had some problems figuring out type errors and precedence rules, but got it pretty soon.

    module Main where
    import System
    
    hailstones :: Integral a => a -> [a]
    hailstones 1 = [1]
    hailstones n
        | odd n = n : hailstones (3*n + 1)
        | even n = n : hailstones (n `div` 2)
    
    main = do
        n <- getArgs
        print . hailstones $ (read  (n !! 0) :: Integer)
    
  2. Graham said

    I recently tried my hand at pulling this off in x86 Assembly; it didn’t go well. Since there’s already a Haskell version, here’s Python:

    #!/usr/bin/env python
    
    def hailstone(n):
        collatz = lambda n: 3 * n + 1 if n & 1 else n // 2
        hs = []
        while n != 1:
            hs.append(n)
            n = collatz(n)
        hs.append(1)
        return hs
    
    if __name__ == "__main__":
        from sys import argv
        print(hailstone(int(argv[1]) if len(argv)> 1 else 1729))
    

    And a more obfuscated C version (posted with the “cpp” sourcecode tag, hope it works):

    #include <stdio.h>
    #include <stdlib.h>
    
    int main(int argc, char **argv) {
        unsigned long n = argc > 1 ? atoi(argv[1]) : 1729;
        while (n != 1) printf("%lu\n", n = (n & 1) ? (n + (n << 1) + 1) : n >> 1);
    }
    
  3. Paul G. said

    And here’s a Perl flavor to display a Hailstone sequence.

    #!/usr/bin/perl -w
    use strict;
    
    print my $n = $ARGV[0];
    print " " . ($n = ($n&1) ? 3*$n+1 : $n/2) while ($n > 1);
    

    I was rather surprised by how fast this works even for the large N. So, I tried figuring out the longest stop length also

    #!/usr/bin/perl -w
    use strict;
    
    my %longest = (4 => 3);
    for my $n (5 .. $ARGV[0]) {
        my $start = $n;
        my @list = ($n);
        push @list, ($n = ($n&1) ? 3*$n+1 : $n/2) while ( ($n > 1) && !exists($longest{$n}) );
        $longest{$start} = scalar(@list);
        $longest{$start} += $longest{$n}-1 if exists($longest{$start});
    }
    
    my $l = 4;
    foreach (keys %longest) {
        $l = $_ if ($longest{$_} > $longest{$l});
    }
    print "longest stop sequence is $longest{$l} for Hailstone($l)\n";
    
  4. ardnew said

    darn, looks like Graham already beat me to the same solution

    #include <stdio.h>
    #include <inttypes.h>
    
    int main(int argc, char *argv[])
    {
      uintmax_t n;
      
      if (argc > 1 && (n = strtoumax(argv[1], 0, 0)) > 0)
      {
        printf("%ju\n", n);
        while (n >> 1)
          printf("%ju\n", n = n & 1 ? n + (n << 1 | 1) : n >> 1);
          
        return 0;
      }
      else
      {
        printf("\nusage:\n\t%s <START>\n\n", argv[0]);
        
        return 1;
      }
    }
    
  5. Mike said

    My 2 cents:

    def hailstone_seq(n):
        """generator for hailstone sequence starting at n
    
        >>>list(hailstone_seq(13))
        [13, 40, 20, 10, 5, 16, 8, 4, 2, 1]
        """"
        yield n
        while n != 1:
            n = (3*n + 1) if n&1 else (n//2)
            yield n
    
    
  6. def hailstone(x):
    “”” A wondrous number is one which eventually reaches 1 using this function.
    See GEB pages 400 and 401. “””
    seq = [x]
    while x != 1:
    # how I miss tail-call optimization in Python
    if not x % 2:
    x /= 2
    else:
    x = 3 * x + 1
    seq.append(x)
    return seq

    # Test
    print hailstone(15)

    # Output
    # [15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]

  7. (Bad formatting above, damn. Now better)

    def hailstone(x):
    """ A wondrous number is one which eventually reaches 1 using this function.
    See GEB pages 400 and 401. """
    seq = [x]
    while x != 1:
    # how I miss tail-call optimization in Python
    if not x % 2:
    x /= 2
    else:
    x = 3 * x + 1
    seq.append(x)
    return seq

    # Test
    print hailstone(15)

    # Output
    # [15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]

  8. slabounty said

    In ruby …

    def hailstone(n)
        begin
            yield n
            n = n % 2 == 1 ? 3*n+1 : n/2
        end while n != 1
        yield n
    end
    
    a = []
    hailstone(13) { |n| a << n }
    p a
    
  9. Eugene said

    This was a rather simple program to solve in Java.

    package hailstones;
    import java.util.Scanner;
    
    public class Hailstones {
    
        public static void main(String[] args) {
            int num1;  //Input variable.
            int count = 0;
            Scanner scan = new Scanner(System.in);
            
            System.out.println("We are going to be running a sequence!");
            System.out.println("Please enter your first number: ");
            num1 = scan.nextInt();
            
            while(num1 != 1)
            {
                if(num1%2 == 0)
                {
                    num1 = num1/2;
                }
                
                else{
                    num1 = (3*num1)+1;
                }
                
                count++;
                System.out.println(num1);
            }
            
            System.out.println("There were "+count+" elements in the sequence.");
        }
    }
    
    
  10. Mike said

    Currently learning racket:

    
    (define (seq n)
      (match n 
        [1 '(1)]
        [_ (cons n (if (even? n) 
                        (seq (/ n 2))
                        (seq (+ (* 3 n) 1))))]))
    
    
  11. Joe Taylor said

    I wass interested in developing a count for the number of items in a sequence for some set of numbers, e.g.:
    1 => 1, length 1
    2 => 2, 1, length 2
    3 => 3, 16, 8, 4, 2, 1, length = 6
    4 => 4, 2, 1, length = 3
    etc.

    The problem seemed to scream “recursion”, and so it can be solved… up to a point, depending on your resources: this code runs in roughly 400ms on a particularly beefy 64-bit desktop with 16GB RAM and quad core, .Net 4.0
    I’m going to see about being a bit less recursive but still performant:
    using System;
    using System.Collections.Generic;
    using System.Diagnostics;
    using System.Linq;
    using System.Text;

    namespace Hailstone
    {
    class Program
    {
    static void Main(string[] args)
    {
    var maxValue = 100000;
    IDictionary stones = new SortedDictionary();
    var stopwatch = new Stopwatch();
    stopwatch.Start();
    //Enumerable.Range(1, maxValue).AsParallel().ForAll(i =>
    // {
    // _stones[i] = Hail(i);
    // });

    foreach (var i in Enumerable.Range(1, maxValue))
    {
    if (!stones.ContainsKey(i))
    {
    stones[i] = Hail(i, stones);
    }
    }

    var byValue = stones.OrderBy(kv => kv.Value).ToDictionary(kv => kv.Key, kv=>kv.Value);
    stopwatch.Stop();
    Console.WriteLine(stopwatch.Elapsed);
    }

    static private int Hail(int x, IDictionary stones)
    {
    if (x == 1)
    {
    stones[x] = x;
    }
    else if (!stones.ContainsKey(x))
    {
    stones[x] = (x % 2 == 0) ?
    1 + Hail(x >> 1, stones) :
    1 + Hail(x*3 + 1, stones);
    }
    return stones[x];
    }
    }
    }

  12. Alec Gorge said
    var s = require('util').puts, n = 13;
    for(;n!=1;n=n&1?3*n+1:n/2,s(n+" "));
    
  13. Joe Taylor said

    OK, so now I have two versions, the recursive, which will very quickly lead to resource abuse, and a non-recursive, which does a little bit less so. Was able to determine paths for fourch 2.1 * 10^6 elements in 4.5 seconds on my seemingly beefy desktop.
    [Apologies for the lack of code style]

    static private int Hail(int x, IDictionary stones)
    {
    int result = 0;
    var trail = new Stack();
    while(x > 0)
    {
    if (stones.ContainsKey(x))
    {
    result = stones[x];
    break;
    }
    else if (x > 1)
    {
    trail.Push(x);
    x = (x%2) == 0 ? x >> 1 : x*3 + 1;
    }
    else // x = 1
    {
    result = stones[x] = 1;
    break;
    }
    }

    while(trail.Count > 0)
    {
    var key = trail.Pop();
    stones[key] = ++result;
    }

    return result;

    }

    static private int RecursiveHail(int x, IDictionary stones)
    {
    if (x == 1)
    {
    stones[x] = 1;
    }
    else if (!stones.ContainsKey(x))
    {
    stones[x] = (x % 2 == 0) ?
    1 + Hail(x >> 1, stones) :
    1 + Hail(x * 3 + 1, stones);
    }
    return stones[x];
    }

  14. Joe Taylor said


    // just checking
    Console.WriteLine("Hello world);

  15. Joe Taylor said

    Caveat:
    Oh my, the path length from 2,139,935,758 is apparently 2 — sometimes 32 bits (unsigned) isn’t enough.

  16. Shankar Chakkere said

    In Apple I BASIC on Mac OS X (http://www.pagetable.com/?p=35)

    #!/usr/bin/apple1basic
    1 REM February 21, 2012 11:25 PM
    100 PRINT “ENTER THE NUMBER”
    200 INPUT N
    300 PRINT N
    310 IF N = 1 THEN END
    410 IF (N/2)*2 = N THEN GOTO 600
    500 REM EVEN NUMBER
    510 N = 3 * N +1
    530 GOTO 300
    600 REM ODD NUMBER
    610 N = N /2
    630 GOTO 300
    900 END

    new-host-4:Apple BASIC shankara_c$ apple1basic hailstone.bas
    ENTER THE NUMBER
    ?13
    13
    40
    20
    10
    5
    16
    8
    4
    2
    1

  17. Manoj Senguttuvan said

    PHP CLI:

    <?
    $xi=$argv[1];
    while($xi!=1)
            echo " ".$xi=$xi%2==0?$xi/2:3*$xi+1;
    ?>
    
  18. This problem lends itself to a recursive solution. Here’s a Javascript snippet exemplifying this fact:

    var hailstone = function(n){      
        if(n === 1)
            return "";    
        
        var cycle = "";
        
        if(n % 2 === 0){
            cycle += (n / 2) + " " + hailstone(n / 2);
        }
        else{
            cycle += (3 * n + 1) + " " + hailstone(3 * n + 1);
        }
        
        return cycle;
    };
    
    console.log(hailstone(13));
    console.log(hailstone(25));
    console.log(hailstone(600));
    

    And here’s the output:

    40 20 10 5 16 8 4 2 1 
    76 38 19 58 29 88 44 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1 
    300 150 75 226 113 340 170 85 256 128 64 32 16 8 4 2 1 
    

    I’m not using memoization or any other cache mechanism, which would allow us to compute larger sequences faster, because I feel that mars the simplicity of the code. Of course, if speed and space were of any real consequence, then that’s an approach that’d have to be explored.

  19. Supriya said

    def hailstone_seq(x):
    print x,
    while x > 1:
    if x % 2 == 0:
    x = x/2
    else:
    x = 3*x + 1
    print x,

  20. Sam Kennedy said

    This is what I made, I want to try and make a fractal but don’t know enough about complex numbers.. I’ll see what I end up with :)

    using namespace std;
    
    int Collatz(int start);
    
    int main(){
    	int i;
    	cout << "Enter Starting Number: ";
    	cin >> i;
    	Collatz(i);
    	system("pause");
    	return 0;
    }
    
    int Collatz(int start){
    	cout << start << " ";
    	if(start == 1){
    		cout << endl;
    		return 1;
    	}
    	if(start % 2 == 0){
    		Collatz(start/2);
    	}
    	else
    	{
    		Collatz((3*start)+1);
    	}
    }
    

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 635 other followers

%d bloggers like this: