Meertens数 in C++ DouKaku?

Meertens数

お題#100「正整数のゲーデル数化？」で定義した goedel を適用すると自分自身になるような数，すなわち goedel (n) == n となるような正整数 n を見つける関数を定義してください．

このような数のことをMeertens数と言うそうです．

32bitsの符号なし整数(あるいは10進10桁整数)までの範囲で探すのにどのくらい計算時間がかかったかをCPUのスペックとともに教えてください．また，その実装で64bit符号なし整数(あるいは10進20桁整数)までの範囲で探すのにどのくらい計算時間がかかりそうか見積ってください(もちろん実際に計算して計算時間を示していただくのでもかまいません)．

補足：

自分でやってみたHaskellでは30行ほどのプログラムで、10進10桁の範囲までの探索は ghci で 4.56秒(CPU Intel Core 2 Duo 2.6GHz)でした。

post comment

Posted feedbacks - C++

多倍長演算を意識する・しないかで言語間の格差が出ますなぁ。

各桁の素数を0～9で乗算したテーブル作成＋論理演算＋枝切り。残念ながら20桁だとunsigned long long(64bit)を溢れてしまうので、19桁までしか計算できません。これ以上は面倒なので勘弁してくださいorz

10桁の計算結果は、1つ(81312000)。計算時間は23.203[秒]（BITSET_SIZE=136, Pen4 3.06GHz, WinXP SP2, VS2005）。ゲーデル数チェック時に使用するbitset(のbit数)に依存する割合が高いので、適正なbit数にすればちょびっと速くなります（10桁を72bitで計算して10.156[秒]）。

枝切りもあるので断定できませんが、テーブル生成時間がほぼ0なので 20桁は10桁の時間×（10の10乗）ぐらいなのかなと。

#include <time.h>
#include <iostream>
#include <bitset>

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

#define BITSET_SIZE (64*2+8)
//#define BITSET_SIZE (32*2+8)

static unsigned long long digit_base[] = {
0x0000000000000001, 0x000000000000000a,
0x0000000000000064, 0x00000000000003E8,
0x0000000000002710, 0x00000000000186A0,
0x00000000000F4240, 0x0000000000989680,
0x0000000005F5E100, 0x000000003B9ACA00,
0x00000002540BE400, 0x000000174876E800,
0x000000E8D4A51000, 0x000009184E72A000,
0x00005AF3107A4000, 0x00038D7EA4C68000,
0x002386F26FC10000, 0x016345785D8A0000,
0x0DE0B6B3A7640000, 0x8AC7230489E80000 };

static unsigned long long prime[20] = {
 2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71 };

static std::bitset<BITSET_SIZE> PRIME_TABLE[ARRAYSIZE(prime)][10];

int CompareBitset(std::bitset<BITSET_SIZE> a, std::bitset<BITSET_SIZE> b) {
    for (int n = BITSET_SIZE - 1; n >= 0; n--) {
        if (a[n] != b[n]) { return a[n] ? 1 : -1; }
    }
    return 0;
}

std::bitset<BITSET_SIZE> bitset_add(std::bitset<BITSET_SIZE> a, std::bitset<BITSET_SIZE> b) {
    std::bitset<BITSET_SIZE> c = (a & b) << 1;
    if (c.any()) { return bitset_add((a ^ b), c); }
    else         { return (a ^ b); }
}

std::bitset<BITSET_SIZE> bitset_multi(std::bitset<BITSET_SIZE> a, std::bitset<BITSET_SIZE> b) {
    std::bitset<BITSET_SIZE> ans;
    if (a.any() && b.any()) {
        for (int i = 0; a.any() && (i < BITSET_SIZE); i++, a <<= 1) {
            if (b[i]) { ans = bitset_add(ans, a); }
        }
    }
    return ans;
}

std::bitset<BITSET_SIZE> bitset_pow(std::bitset<BITSET_SIZE> a, unsigned int n) {
    std::bitset<BITSET_SIZE> v((unsigned long)1);
    while (n != 0) {
        if (n & 0x1)       { v = bitset_multi(v, a); }
        if ((n /= 2) != 0) { a = bitset_multi(a, a); }
    }
    return v;
}

unsigned int digit(unsigned long long n, unsigned int base)
{
    unsigned int d = 1;
    for (; n > (base - 1); n /= base) { d++; }
    return d;
}

std::bitset<BITSET_SIZE> ulonglong2bitset(unsigned long long v) {
    std::bitset<BITSET_SIZE> high((unsigned long)(v >> (sizeof(v)*4)));
    std::bitset<BITSET_SIZE> low((unsigned long)v);
    return ((high << (sizeof(v)*4)) | low);
}

void InitPrimeTable() {
    for (int i = 0; i < ARRAYSIZE(PRIME_TABLE); i++) {
        for (int j = 0; j < ARRAYSIZE(PRIME_TABLE[0]); j++) {
            PRIME_TABLE[i][j] = bitset_pow(ulonglong2bitset(prime[i]), j);
        }
    }
}

int CheckMeertensNumber(unsigned long long n) {
    std::bitset<BITSET_SIZE> m = ulonglong2bitset(n);
    std::bitset<BITSET_SIZE> g(1UL);
    for (unsigned int d = digit(n, 10); d > 0; d--, n /= 10) {
        unsigned int p = (unsigned int)(n % 10);
        if (CompareBitset(PRIME_TABLE[d - 1][p], m) > 0)
            return 1;
        g = bitset_multi(g, PRIME_TABLE[d - 1][p]);
        if (CompareBitset(g, m) > 0)
            return 1;
    }
    return (m == g) ? 0 : -1;
}

int SearchMeertensNumber(unsigned long long base, unsigned int d, int& cnt) {
    if (d <= 0) {
        int ret = CheckMeertensNumber(base);
        if (ret == 0) {
            std::cout << base << std::endl;
            cnt++;
        }
        return ret;
    }
    if (SearchMeertensNumber(base, d - 1, cnt) > 0)
        return (base == 0) ? 0 : 1;
    for (unsigned long long i = 1; i < 10; i++) {
        base += digit_base[d - 1];
        if (SearchMeertensNumber(base, d - 1, cnt) > 0)
            break;
    }
    return 0;
}

void main(int argc, char* argv[]) {
    unsigned int nDigit = 10;
    int nCount = 0;

    std::cout << "---" << std::endl;
    std::cout << "bitsize=" << BITSET_SIZE << std::endl;
    std::cout << "digit  =" << nDigit << std::endl;

    clock_t st = clock();
    InitPrimeTable();
    clock_t mt = clock();
    SearchMeertensNumber(0, nDigit, nCount);
    clock_t et = clock();

    std::cout << "count=" << nCount << std::endl;
    std::cout << "init[s]:" << ((double)(mt-st)/CLOCKS_PER_SEC) << std::endl;
    std::cout << "calc[s]:" << ((double)(et-mt)/CLOCKS_PER_SEC) << std::endl;
    std::cout << "full[s]:" << ((double)(et-st)/CLOCKS_PER_SEC) << std::endl;
}

コンパイラの最適化具合によるかもしれませんが、リファクタリングによって少し速度が改善されたので訂正。

bitサイズ136で、10桁 5.25[秒]，11桁 19.75[秒]，12桁 66.343[秒]，13桁 222.078[秒]。（ちなみにbitサイズ72で、10桁 3.469[秒]）

どうも1桁あたり10倍ではなく、3.5～4.0倍になってるようなので、20桁で18日～60日程度と予想（幅ありすぎ）。

//ソースは置換する関数のみ
#include <limits>

int CompareBitset(std::bitset<BITSET_SIZE> a, std::bitset<BITSET_SIZE> b) {
    std::bitset<BITSET_SIZE> mask( std::numeric_limits<unsigned long>::max() );
    int ans = 0;
    while (a.any() || b.any()) {
        unsigned long an = (a & mask).to_ulong();
        unsigned long bn = (b & mask).to_ulong();
        if      (an > bn) { ans =  1; }
        else if (an < bn) { ans = -1; }
        a >>= std::numeric_limits<unsigned long>::digits;
        b >>= std::numeric_limits<unsigned long>::digits;
    }
    return ans;
}

std::bitset<BITSET_SIZE> bitset_multi(std::bitset<BITSET_SIZE> a, std::bitset<BITSET_SIZE> b) {
    std::bitset<BITSET_SIZE> ans;
    for (;a.any() && b.any(); b >>= 1, a <<= 1) {
        if (b[0]) { ans = bitset_add(ans, a); }
    }
    return ans;
}

DouKaku?

Meertens数

Posted feedbacks - C++

Index

Feed

Other

Link

Pathtraq