Project Euler 850

题目

Fractions of Powers

Any positive real number \(x\) can be decomposed into integer and fractional parts \(\lfloor x \rfloor + \{x\}\), where \(\lfloor x \rfloor\) (the floor function) is an integer, and \(0\le \{x\} < 1\).

For positive integers \(k\) and \(n\), define the function

\[\begin{aligned} f_k(n) = \sum_{i=1}^{n}\left\{ \frac{i^k}{n} \right\} \end{aligned}\]

For example, \(f_5(10)=4.5\) and \(f_7(1234)=616.5\).

Let

\[\begin{aligned} S(N) = \sum_{\substack{k=1 \\k\text{ odd}}}^{N} \sum_{n=1}^{N} f_k(n) \end{aligned}\]

You are given that \(S(10)=100.5\) and \(S(10^3)=123687804\).

Find \(\lfloor S(33557799775533) \rfloor\). Give your answer modulo \(977676779\).

解决方案

固定 \(n\)，令 \(k\) 为奇数。对 \(1\le i\le n-1\)，考虑配对 \((i,n-i)\)。因为 \(k\) 是奇数，\((n-i)^k\equiv (-i)^k\equiv -i^k\pmod n.\)令 \(i^k=qn+r\)，其中 \(0\le r<n\)。则\(\left\{\dfrac{i^k}{n}\right\}=\dfrac{r}{n}.\)而 \((n-i)^k\equiv -r\pmod n\)，即其余数为 \(0\) 当且仅当 \(r=0\)，否则余数为 \(n-r\)，于是

\[ \left\{\frac{(n-i)^k}{n}\right\}= \begin{cases} 0,& r=0,\\ \frac{n-r}{n},& r\ne 0. \end{cases} \]

因此对 \(1\le i\le n-1\)， \[ \left\{\frac{i^k}{n}\right\}+\left\{\frac{(n-i)^k}{n}\right\}=\mathbf{1}\{n\nmid i^k\} \]

另外当 \(i=n\) 时 \(\{n^k/n\}=0\)。

令\(g_k(n)=\#\{1\le i\le n: n\mid i^k\}.\)则在 \(i=1,\dots,n\) 中，满足 \(n\mid i^k\) 的项贡献 \(0\)，其余项在配对后总贡献为 \(1\)，从而\(f_k(n)=\dfrac{n-g_k(n)}{2} (k\equiv 1\pmod 2).\)于是

\[ S(N)=\sum_{\substack{k=1\\k\equiv 1\pmod 2}}^{N}\sum_{n=1}^{N}f_k(n) =\frac12\left(\Big\lfloor\frac{N+1}{2}\Big\rfloor\sum_{n=1}^{N}n-\sum_{\substack{k=1\\k\equiv 1\pmod 2}}^{N}G_k(N)\right), \]

其中\(\displaystyle G_k(N)=\sum_{n=1}^{N}g_k(n).\)

注意到 \(2S(N)\) 是整数（因为 \(2f_k(n)=n-g_k(n)\) 为整数），因此计算 \(\lfloor S(N)\rfloor\bmod M\) 时，可先算\(A=2S(N)\bmod{2M},\)再由 \[ \lfloor S(N)\rfloor= \begin{cases} A/2,&A\equiv 0\pmod 2,\\ (A-1)/2,&A\equiv 1\pmod 2 \end{cases} \] 在模 \(M\) 下恢复答案（因为模数取为 \(2M\)，整除 \(2\) 在整数意义下可直接做）。

若 \(\gcd(a,b)=1\)，由中国剩余定理，\(i\bmod ab\) 与 \((i\bmod a,i\bmod b)\) 一一对应，并且\(ab\mid i^k \iff a\mid i^k\land b\mid i^k.\)因此解的个数可以相乘，\(g_k\) 为积性函数。

对 \(n=p^e\)，条件 \(p^e\mid i^k\) 等价于 \(v_p(i^k)\ge e\)，即\(kv_p(i)\ge e\iff v_p(i)\ge \left\lceil\dfrac{e}{k}\right\rceil.\)在 \(1\le i\le p^e\) 中，满足 \(v_p(i)\ge t\) 的数正是 \(p^t\) 的倍数，共 \(p^{e-t}\) 个。取 \(t=\lceil e/k\rceil\) 得\(g_k(p^e)=p^{e-\lceil e/k\rceil} (e\ge 1).\)

令常值函数 \(\mathbf{1}(n)=1\)，其Dirichlet逆为莫比乌斯函数 \(\mu\)，即 \(\mathbf{1}^{-1}=\mu\)。定义\(h_k = g_k \ast \mu,\)则\(g_k = \mathbf{1} \ast h_k.\) 对任意 \(N\)， \[ G_k(N)=\sum_{n\le N}g_k(n) =\sum_{n\le N}\sum_{d\mid n}h_k(d) =\sum_{d\le N}h_k(d)\Big\lfloor\frac{N}{d}\Big\rfloor. \]

由 \(h_k=g_k*\mu\) 的乘法性可得 \(h_k\) 也是积性函数，并且对素幂\(h_k(p)=g_k(p)-g_k(1)=1-1=0,h_k(p^e)=g_k(p^e)-g_k(p^{e-1}) (e\ge 2).\)因此当且仅当每个质因子的指数都不小于 \(2\) 时，\(h_k(n)\) 才可能非零；这类整数称为 powerful 数。

于是计算 \(G_k(N)\) 只需枚举所有 \(d\le N\) 的powerful 数并累加 \(h_k(d)\lfloor N/d\rfloor\)。powerful 数个数为 \(O(\sqrt N)\)，可用 DFS 按质因子递增枚举得到。

对任意 \(n\le N\)，若 \(p^e\mid n\)，则 \(p^e\le N\)，从而\(e\le \log_p N \le \log_2 N.\)若 \(k>\log_2 N\)，则对所有出现的 \(e\) 都有 \(e<k\)，于是 \(\lceil e/k\rceil=1\)，从而\(g_k(p^e)=p^{e-1},\)完全与 \(k\) 无关，故 \(G_k(N)\) 在所有足够大的 \(k\) 上稳定为常数。因为这里只对奇数 \(k\) 求和，只需取\(k_0=\min\{k:k\equiv1\pmod 2\land k>\lfloor \log_2 N\rfloor\},\)则对所有奇数 \(k\ge k_0\) 有 \(G_k(N)=G_{k_0}(N)\)。因此

\[ \sum_{\substack{k=1\\k\equiv1\pmod 2}}^{N}G_k(N) =\sum_{\substack{k<k_0\\k\equiv1\pmod 2}}G_k(N) +\left(\Big\lfloor\frac{N+1}{2}\Big\rfloor-\frac{k_0-1}{2}\right)G_{k_0}(N). \]

最终计算\(A=2S(N)\)，回代计算出\(S(N)\)的最终值即可。

代码

#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;
const ll N = 33557799775533LL;
const ll M = 977676779LL;
const ll MOD2 = 2LL * M;

vector<int> primes;
int cur_k;

static inline ll mul2(ll a, ll b)
{
    return (ll)((a % MOD2) * (b % MOD2) % MOD2);
}

ll dfs(int start, ll ncur, ll vcur, ll N0)
{
    ll res = mul2(vcur, (N0 / ncur) % MOD2);
    ll limit = N0 / ncur;

    for (int idx = start; idx < (int)primes.size(); ++idx)
    {
        ll p = primes[idx];
        if (p * p > limit)
            break;

        ll powp = p * p;
        int e = 2;
        ll gcur = 1 % MOD2;

        while (powp <= limit)
        {
            ll gprev = gcur;
            if (e % cur_k != 1)
                gcur = mul2(gcur, p % MOD2);
            ll hk = (gcur - gprev + MOD2) % MOD2;

            ll nnext = ncur * powp;
            ll vnext = mul2(vcur, hk);

            if (vnext)
            {
                res += dfs(idx + 1, nnext, vnext, N0);
                res %= MOD2;
            }

            if (powp > limit / p)
                break;
            powp *= p;
            ++e;
        }
    }

    return res % MOD2;
}

ll Gk(int k, ll N0)
{
    if (k == 1)
        return N0 % MOD2;
    cur_k = k;
    return dfs(0, 1LL, 1LL, N0);
}

ll solve(ll N0)
{
    ll lim = int_square(N0);
    primes = Factorizer(lim).primes;

    int emax = 63 - __builtin_clzll((unsigned long long)N0);
    int k0 = (emax % 2 == 0) ? (emax + 1) : (emax + 2);

    ll K = (N0 + 1) / 2;

    ll sum_n_mod2;
    if ((N0 & 1LL) == 0)
        sum_n_mod2 = mul2((N0 / 2) % MOD2, (N0 + 1) % MOD2);
    else
        sum_n_mod2 = mul2(((N0 + 1) / 2) % MOD2, (N0 % MOD2));

    ll term1 = mul2(K % MOD2, sum_n_mod2);

    ll sumG = 0;

    if ((ll)k0 <= N0)
    {
        for (int k = 1; k <= k0 - 2; k += 2)
        {
            sumG += Gk(k, N0);
            sumG %= MOD2;
        }
        ll gst = Gk(k0, N0);
        ll cnt = K - (ll)((k0 - 1) / 2);
        sumG += mul2(cnt % MOD2, gst);
        sumG %= MOD2;
    }
    else
    {
        for (ll k = 1; k <= N0; k += 2)
        {
            sumG += Gk((int)k, N0);
            sumG %= MOD2;
        }
    }

    ll A = (term1 - sumG + MOD2) % MOD2;

    ll ans;
    if ((A & 1LL) == 0)
        ans = (A / 2) % M;
    else
        ans = ((A - 1) / 2) % M;

    return ans;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    cout << solve(N) << "\n";
    return 0;
}