[Project Euler 521]Smallest prime factor 题解

为什么博客园的LaTex有问题啊。

 跑了40s左右，跑得还是有点久……

//waz
#include <bits/stdc++.h>

using namespace std;

#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define ALL(x) (x).begin(), (x).end()
#define SZ(x) ((int)((x).size()))

typedef pair<int, int> PII;
typedef vector<int> VI;
typedef long long int64;
typedef unsigned int uint;
typedef unsigned long long uint64;

#define gi(x) ((x) = F())
#define gii(x, y) (gi(x), gi(y))
#define giii(x, y, z) (gii(x, y), gi(z))

int F()
{
    char ch;
    int x, a;
    while (ch = getchar(), (ch < '0' || ch > '9') && ch != '-');
    if (ch == '-') ch = getchar(), a = -1;
    else a = 1;
    x = ch - '0';
    while (ch = getchar(), ch >= '0' && ch <= '9')
        x = (x << 1) + (x << 3) + ch - '0';
    return a * x;
}

const int mod = 1e9;

int inc(int a, int b) { a += b; return a >= mod ? a - mod : a; }

int dec(int a, int b) { a -= b; return a < 0 ? a + mod : a; }

const int N = 1e6 + 10;

int64 n;

int sump[N], sum2[N];

int sqrtn;

int p[N], m;

namespace work1
{
    int sum(int64 x)
    {
        int64 y = x + 1;
        if (x & 1) y >>= 1; else x >>= 1;
        return (x % mod) * (y % mod) % mod;
    }
    int S(int64 x, int j)
    {
        if (x < sqrtn && x < 1LL * p[j + 1] * p[j + 1]) return sum2[x];
        if (!j) { return dec(sum(x), 1); }
        if (1LL * p[j] * p[j] > x) return S(x, j - 1);
        else return dec(S(x, j - 1), 1LL * p[j] * dec(S(x / p[j], j - 1), sump[j - 1]) % mod);
    }
}

namespace work2
{
    const int K = 5000;
    unordered_map<int64, int> ha;
    bitset<200010000> s;
    int dfs(int64 x, int i)
    {
        if (!i) return x % mod;
        if (!x) return 0;
        if (ha[i * n + x]) return ha[i * n + x];
        int ans = x % mod;
        for (int j = i; j; --j)
        {
            ans = dec(ans, dfs(x / p[j], j - 1));
        }
        return ha[i * n + x] = ans;
    }
    int solve()
    {
        int ans = 0;
        int k = 1;
        for (; p[k] < K && k < m; ++k);
        //cerr << k << endl;
        for (int i = 1; i <= k; ++i)
            ans = inc(ans, 1LL * p[i] * (dfs(n / p[i], i - 1) - 1) % mod);
        if (k == m) return ans;
        int64 ret = n / p[k] - 1;
        for (int i = 1; i <= k; ++i)
            for (int j = p[i], c = n / p[k]; j <= c; j += p[i])
                if (!s[j]) s[j] = 1, --ret;
        for (int i = k + 1; i <= m; ++i)
        {
            int t = n / p[i - 1], tt = n / p[i];
            for (int j = t; j > tt; --j)
                if (!s[j]) --ret;
            ans = inc(ans, ret % mod * p[i] % mod);
            for (int j = p[i]; j <= tt; j += p[i])
                if (!s[j]) s[j] = 1, --ret;
        }
        return ans;
    }
}

main()
{
    cin >> n;
    sqrtn = sqrt(n);
    static bool fg[N];
    for (int i = 2; i <= sqrtn; ++i)
    {
        if (!fg[i])
            p[++m] = i;
        for (int j = 1; j <= m && p[j] * i <= sqrtn; ++j)
        {
            fg[p[j] * i] = 1;
            if (i % p[j] == 0) break;
        }
    }
    for (int i = 1; i <= m; ++i) sump[i] = inc(sump[i - 1], p[i]);
    for (int i = 2; i <= sqrtn; ++i) sum2[i] = inc(sum2[i - 1], i * (!fg[i]));
    int ans2 = work2::solve();
    //cerr << "task2 : " << clock() << " ms" << endl;
    //cerr << ans2 << endl;
    int ans1 = work1::S(n, m);
    //cerr << "task1 : " << clock() << " ms" << endl;
    //cerr << ans1 << endl;
    cout << inc(ans1, ans2) << endl;
}