Project Euler - Problem 193

Squarefree Numbers

Problem 193 (日本語訳)

N = 2⁵⁰ 以下の無平方数はいくつか？

すぐ思い付くのはエラトステネスの篩っぽく平方数の倍数をふるい落とす方法ですが，N が 2⁵⁰ とでかいので無理です．

そこで包除原理を使う．

無平方数でない(平方因子をもつ) N 以下の正整数の総数を S，
n² で割り切れる N 以下の正整数の個数を S(n) とすると，

S = S(2) + S(3) + S(5) + … - S(2･3) - S(2･5) - S(3･5) - … + S(2･3･5) + … - S(2･3･5･7) - …

これを sqrt(N) = 2²⁵ までの素数にわたり計算します．組合せの総数がヤバいことになりそうですが

2･3･5･7･11･13･17･19･23 < 2²⁵ < 2･3･5･7･11･13･17･19･23･29

で積は高々 9 個までなのでそこまでヤバいことにはなりません．

typedef long long ll;

ll N = (ll)1 << 25;     // 2^25
vector<ll> primes;        // 2^25 以下の素数
ll S = 0;              // 無平方数でない数の総数

// prod : 掛け合わせてきた素数の積
// pos : 次に掛ける素数の位置
// sgn : 奇数個の積なら 1, 偶数個なら -1
void rec(ll prod, int pos, int sgn)
{
    ll n;
    // 積が 2^25 = sqrt(2^50) を超えない範囲で
    while (pos < primes.size() && (n = prod * primes[pos]) <= N) {
        // S(n) = floor(2^50 / n^2)
        // 奇数個の積なら加え，偶数個の積なら引く
        S += N * N / (n * n) * sgn;
        // 素数集合のポインタをインクリメント，符号を反転して再帰
        rec(n, ++pos, -sgn);
    }
}

int PE0193()
{
    // ふるい
    bool* p = new bool[N + 1];
    p[0] = p[1] = false;
    for (int i = 2; i <= N; i++) p[i] = true;

    for (int i = 2; i * i <= N; i++)
        if (p[i])
            for (int j = i * i; j <= N; j += i)
                p[j] = false;

    // primes : 2^25 以下の素数の集合
    for (int i = 0; i <= N; i++)
        if (p[i])
            primes.push_back(i);

    // 再帰により S を計算
    rec(1, 0, 1);

    // answer
    cout << (N * N) - S << endl;

    return 0;
}


int main()
{
    clock_t start, end;
    start = clock();

    PE0193();

    end = clock();
    printf("%d msec.\n", end - start);

    return 0;
}

20 秒ちょいかかった．