DouKaku? 4872 ところてん: とりあえず、再帰で書いてみた。 sea...(Meertens数)

Comment detail

とりあえず、再帰で書いてみた。
searchMeertens(n)でn桁の範囲で探索します。
8桁で一つ見つかりました。

8桁を調べるのに160秒かかりました。
桁数が上がるごとに10倍の計算量になるので、
10桁では4時間くらいかなぁ。

とりあえず、累乗計算をしないようにしたけど、これからどうしようかな。
もうちょっと枝が刈れそうな気がするけど。

Cで書けばもうちょっと早くなるかな。

import time

def CreatePrimes(n):
    primes = []
    c = 2
    while len(primes) < n:
        primes.append(c)
        for x in primes[:-1]:
            if c%x == 0:
                primes.pop();
                break
        c += 1
    return primes

def searchMeertens(keta):
    def _search(keta,count,r,n,primes,first):
        n *= 10
        for x in range(first,10):
            if x:
                r *= primes[count]
                n += 1
            if r == n:
                print "Found Meertens.",n
            if keta > 1:
                _search(keta - 1, count + 1, r ,n ,primes, 0)
    p = CreatePrimes(20)
    _search(keta,0,1,0,p,1)


keta = 10
t0 = time.time()
t1 = t0

for n in range(1,keta+1):
    searchMeertens(n)
    print n , time.time() - t1 , "sec"
    t1 = time.time()
    
print "total time", time.time() - t0

Cで書き直したら恐ろしく速くなった。 10桁実行するのに176秒＠AthlonXP2500+ pythonで書いたのと比較して、だいたい100倍高速化した計算に。

ゲーテル数って10桁でintの範囲をあっさりと突破するのね・・・バグに気づかなかった。 unsigned long long int なんて初めて使ったよ。

#include <stdio.h>
#include <time.h>

int primes[] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};

void searchMeertens(int keta, int count, unsigned long long int r, unsigned long long int n, int first){
    int x;
    n *= 10;
    for(x = first; x < 10; x++){
        if (x){
            r *= primes[count];
            n++;
        }    
        if (r == n){
            printf("Found Meertens %d\n",n);
        }
        if (keta > 1){
            searchMeertens(keta -1, count +1, r, n, 0);
        }
    }
}


int main(void)
{
    int keta;
    int x;
    time_t t;
    keta = 10;
    t = time(NULL);
    for(x = 1; x <= keta; x++){
        searchMeertens(x,0,1,0,1);
        printf("%d %dsec\n",x, time(NULL)-t);
        t = time(NULL);
    }
    printf("push anykey to exit.");
    getchar();
}

一応手元で11桁まで回り終わりましたが、
一桁増えると計算時間も10倍になりました。
(当然と言えば当然ですが)

なので、google先生に聞いたところ、20桁では
(176 * (10^10)) *秒 = 55 772.2258 年
だそうです。

とりあえず40秒くらいまで縮まったけど、
他の人が何で爆速なのかが分からない。
あとどんな枝狩りすりゃいいねん。

枝狩りで分岐予測がこける関係で、Pen4系では遅いと思う。

関数型言語は読めないよ・・・

#include <stdio.h>
#include <time.h>

int primes[] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};

void _searchMeertens(int keta, int count, unsigned long long int r, unsigned long long int n){
    int x;
    n *= 10;
    for(x = (count == 0); x < 10; x++){
        if (x){
            r *= primes[count];
            n++;
        }    
        if (keta == 0){
            if(r > n){        //枝狩り
                return;    
            }
            else if (r == n){
                printf("Found Meertens %d\n",n);
            }
        }else{
            _searchMeertens(keta -1, count +1, r, n);
        }
    }
}

void searchMeertens(int keta)
{
    int x;
    time_t t;
    t = time(NULL);
    for(x = 0; x < keta; x++){
        _searchMeertens(x,0,1,0);
        printf("keta = %d time = %d\n", x+1, time(NULL) - t);
        t = time(NULL);
    }
}

int main(void)
{
    int keta;
    time_t t;

    keta = 10;
    t = time(NULL);
    searchMeertens(keta);
    printf("%dsec\n", time(NULL)-t);

    printf("push anykey to exit.");
    getchar();
}

r が、その桁の最大値を超えたときを狩れば結構削れますよ。

void _searchMeertens(int keta, int count, unsigned long long int r, unsigned long long int n, unsigned long long int nmax)
{
    int x;
    n *= 10;
    for (x = (count == 0); x < 10; x++) {
        if (x) {
            r *= primes[count];
            if (r > nmax)    /* 枝狩り2 */
                return;
            n++;
        }
        if (keta == 0) {
            if(r > n)       /* 枝狩り*/
                return;
            if (r == n)
                printf("Found Meertens %d\n", n);
        }else{
            _searchMeertens(keta -1, count + 1, r, n, nmax);
        }
    }
}

void searchMeertens(int keta)
{
    int x;
    unsigned long long int nmax;
    time_t t;
    t = time(NULL);
    for (x = 0; x < keta; x++) {
        nmax = (unsigned long long)pow(10.0, x + 1) - 1UL;
        _searchMeertens(x, 0, 1, 0, nmax);
        printf("keta = %d time = %d\n", x + 1, time(NULL) - t);
        t = time(NULL);
    }
}

なるほどなーと納得しました。

手元の奴で10桁が1秒未満で完了した・・・。
確かにそうすればよかったのかと後悔してます。

ちなみに私の環境での実行結果。
あまり時間が無かったので、18桁まで実行させましたが、20桁は20分くらいですね。
なのでトータルタイムは30分くらいでしょうか。

keta = 1 time = 0
keta = 2 time = 0
keta = 3 time = 0
keta = 4 time = 0
keta = 5 time = 0
keta = 6 time = 0
keta = 7 time = 0
Found Meertens 81312000
keta = 8 time = 0
keta = 9 time = 0
keta = 10 time = 0
keta = 11 time = 0
keta = 12 time = 0
keta = 13 time = 1
keta = 14 time = 2
keta = 15 time = 6
keta = 16 time = 16
keta = 17 time = 45
keta = 18 time = 127

最初の乗算(n *= 10)を外に出してみれば、もうちょいいけるかも....。

keta = 14 time = 2
keta = 15 time = 4
keta = 16 time = 12
keta = 17 time = 34
keta = 18 time = 96

微々たるもんしか減らないorz。

#include <stdio.h>
#include <time.h>

#ifndef ARRAYSIZE
#define ARRAYSIZE(x)  (sizeof(x)/sizeof(x[0]))
#endif

unsigned long long int dec_base[20];

void _searchMeertens(int keta, int d, unsigned long long int r, unsigned long long int n) {
    int x;
    for (int x = (d == 0); x < 10; x++) {
        if (x > 0) {
            r *= primes[d];
            if (r >= dec_base[keta])
                return;
            n += dec_base[keta - (d + 1)];
        }
        if ((d + 1) >= keta) {
            if (r >= n) {
                if (r == n)
                    printf("Found Meertens %lld\n", n);
                return;
            }
        } else {
            _searchMeertens(keta, d + 1, r, n);
        }
    }
}

void searchMeertens(int keta) {
    int x;
    time_t t;
    t = time(NULL);
    for (x = 1; x <= keta; x++) {
        _searchMeertens(x, 0, 1, 0);
        printf("keta = %d time = %d\n", x, time(NULL) - t);
        t = time(NULL);
    }
}

void main() {
    int i;
    dec_base[0] = 1;
    for (i = 1; i < ARRAYSIZE(dec_base); i++) {
        dec_base[i] = dec_base[i - 1] * 10;
    }
    searchMeertens(19);
}

DouKaku?

Comment detail

Index

Feed

Other

Link

Pathtraq