いつの魔にRound 2がはじまっている。

LLR2006
 キミならどう書く 2.0 - ROUND 2 - ? Lightweight Language Ring
お題は「Collatz予想」（角谷予想，3x+1問題）についての問題です．
[中略]
また，余力のある者は大きな n についても h(n) を求めよ．

これまたtrivialな問題なのだけど、1,000,000とか計算する場合、実は32bit整数では足りない。

タイトルどおり、113383で引っかかるのだ。

113383 → 340150 → 170075 → 510226 → 255113 → 765340 → 382670 → 191335 → 574006 → 287003 → 861010 → 430505 → 1291516 → 645758 → 322879 → 968638 → 484319 → 1452958 → 726479 → 2179438 → 1089719 → 3269158 → 1634579 → 4903738 → 2451869 → 7355608 → 3677804 → 1838902 → 919451 → 2758354 → 1379177 → 4137532 → 2068766 → 1034383 → 3103150 → 1551575 → 4654726 → 2327363 → 6982090 → 3491045 → 10473136 → 5236568 → 2618284 → 1309142 → 654571 → 1963714 → 981857 → 2945572 → 1472786 → 736393 → 2209180 → 1104590 → 552295 → 1656886 → 828443 → 2485330 → 1242665 → 3727996 → 1863998 → 931999 → 2795998 → 1397999 → 4193998 → 2096999 → 6290998 → 3145499 → 9436498 → 4718249 → 14154748 → 7077374 → 3538687 → 10616062 → 5308031 → 15924094 → 7962047 → 23886142 → 11943071 → 35829214 → 17914607 → 53743822 → 26871911 → 80615734 → 40307867 → 120923602 → 60461801 → 181385404 → 90692702 → 45346351 → 136039054 → 68019527 → 204058582 → 102029291 → 306087874 → 153043937 → 459131812 → 229565906 → 114782953 → 344348860 → 172174430 → 86087215 → 258261646 → 129130823 → 387392470 → 193696235 → 581088706 → 290544353 → 871633060 → 435816530 → 217908265 → 653724796 → 326862398 → 163431199 → 490293598 → 245146799 → 735440398 → 367720199 → 1103160598 → 551580299 → 1654740898 → 827370449 → -1812855948

64bitなら、1,000,000まで計算する場合でも何とか足りる。

というわけで、以下、Cによる実装。

#include <stdio.h>
#include <stdlib.h>

typedef long long I64;

#define collatz(n) ((n) == 1 ? 1 : ((n) % 2) ? (n)*3+1 : (n)/2)

I64 g(I64 n){
    I64 i = n;
    I64 result = 1;
    do {
        result++;
    }while((i = collatz(i)) != 1);
    return result;
}

I64 *h(I64 n){
    I64 i, nmax = 0, gmax = 0, gnext = 0;
    static I64 result[2];
    for (i = 1; i <= n; i++){
        gnext = g(i);
#ifdef VERBOSE
        fprintf(stderr, "g(%qd) = %qd\r", i, gnext);
#endif
        if (gnext > gmax){
            nmax = i; gmax = gnext;
        }
    }
    result[0] = nmax; result[1] = gmax;
    return result;
}

int main(int argc, char **argv){
    I64 n = 100;
    if (argc > 1){
        n = atoi(argv[1]);
    }
    I64 *hg = h(n);
    printf("h(%qd) = %qd where g(%qd) = %qd\n", n, hg[0], hg[0], hg[1]);
    return 0; /* so -Wall is happy */
}

MacBook Pro 2GHzでの実行結果は、以下のとおり。

% gcc -Wall -O3 iterative.c
% /usr/bin/time ./a.out 1000000
h(1000000) = 837799 where g(837799) = 525
        3.19 real         3.13 user         0.01 sys

C, C++, Perlは気をつけられたし。Haskell, Ruby, Python は整数はデフォルトでbignumなのでOKこの問題はない。なお、Perlの場合-Duse64bitintを付けてcompileされたものであれば1,000,000も目指せる。use bignumでもOKだけど低速ゆえお薦めしない。

Dan the Integral Man

追記:id:palmoさんのようにビット演算でやる方法もあるが、この場合でもやはり113383で引っかかる。以下、検証コードと結果。

iterative.pl 
use strict;
use warnings;

sub g($) {
    use integer;
    my $num = shift;
    my $n = $num;
    my $step = 1;
    while ($n > 1) {
        $n = ($n & 1) ? (3 * $n + 1) : ($n >> 1);
        $step++;
    }
    die "g($num) got $n!" unless $n == 1;
    return $step;
}

sub h{
    my $n = shift;
    my ($nmax, $gmax, $gnext) = (0, 0, 0);
    for (1..$n){
        my $gnext = g($_);
        print STDERR "g($_) = $gnext\r";
        next unless $gnext > $gmax;
        $nmax = $_;
        $gmax = $gnext;
    }
    return ($nmax, $gmax);
}

sub printh{
    my $n = shift;
    my ($h, $g) = h($n);
    print "h($n) = $h where g($n) = $g\n";
}

my $n = shift || 100;
printh $n;
__END__
実行結果: 
% perl -V | grep 64bit
    use64bitint=undef use64bitall=undef uselongdouble=undef
% perl iterative.pl 1000000
g(113383) got -1812855948! at iterative.pl line 13.

追記^2

をを、itaさんも同じとこでハマってる。

JGeek Log - 3n+1 問題
1000000 でint のオーバーフローで落ちる。