uestc 811 GCD 杜教筛 + 自然幂和

传送门:UESTC 811

---

这题思路出的还算快, 但是UESTC OJ评测卡了一天, 代码细节比较多, 调bug的时候都怀疑解法是不是又错了

这题好像不好找题解

题意

1<=N<=10^10, 1<=K<=5

求$\sum_{i=1}^N\sum_{j=1}^Ngcd(i, j)^k$

---

题解

枚举gcd, 推下式子

$$\begin{align*} ans &= \sum_{i=1}^N\sum_{j=1}^Ngcd(i, j)^k= 2\sum_{i=1}^N\sum_{j=1}^igcd(i, j)^k-\sum_{i=1}^Ngcd(i, i)^k \\ &=2\sum_{i=1}^N\sum_{d|i}d^k\phi(\frac{i}{d})-\sum_{i=1}^Ni^k\\ &=2\sum_{i=1}^Ni^k\sum_{d=1}^{\lfloor\frac{N}{i}\rfloor}\phi(d)-\sum_{i=1}^Ni^k \\ \end{align*}$$

然后就可以分块做了

---

快速求自然幂和

$\sum_{i=1}^Ni^k=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^k$
预处理出伯努利数和组合数就可以在O(k+1)时间内求出自然幂和,
k最大为5

---

伯努利数求法:
$B_0 = 1, B_k=-\frac{1}{k+1}\sum_{i=0}^{k-1}C_{k+1}^iB_i$
求法源于:
$\sum_{i=0}^kC_{k+1}^iB_i=0$

---

欧拉函数的杜教筛之前题解里涉及到, 就不多讲了
总复杂度$O(n^{\frac{2}{3}})$

---

code:

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
#define next Next
const int N = 10000001;
const LL mod = 1e9 + 7;
const LL _2 = (mod + 1) / 2;
const int mo = 2333333;

bool isPrime[N];
int prime[N / 5];
LL phi[N];
LL C[10][10], B[10];

int cnt;

LL qpow(LL a, int b)
{
    a %= mod;
    LL res = 1;
    while(b)
    {
        if(b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

void init()
{
    for(int i = 0; i < 10; ++i)
        C[i][0] = C[i][i] = 1;
    for(int i = 2; i < 10; ++i)
        for(int j = 1; j < i; ++j)
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
    B[0] = 1;
    for(int i = 1; i < 10; ++i)
    {
        if(i == 9) continue;
        LL tmp = 0;
        for(int j = 0; j < i; ++j)
            tmp = (tmp + C[i + 1][j] * B[j]) % mod;
        tmp = tmp * qpow(i + 1, mod - 2) % mod;
        B[i] = (mod - tmp) % mod;
    }
    memset(isPrime, true, sizeof isPrime);
    phi[1] = 1;
    isPrime[1] = false;
    for(int i = 2; i < N; ++i)
    {
        if(isPrime[i])
        {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for(int j = 1; j <= cnt && i * prime[j] < N; ++j)
        {
            isPrime[i * prime[j]] = false;
            if(i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
    for(int i = 2; i < N; ++i)
        phi[i] = (phi[i] + phi[i - 1]) % mod;
}

LL SumOfPow(LL n, int k)
{
    LL res = 0;
    for(int i = 1; i <= k + 1; ++i)
        res = (res + C[k + 1][i] * B[k + 1 - i] % mod * qpow(n + 1, i) % mod) % mod;

    LL inv = qpow(k + 1, mod - 2);
    res = res * inv % mod;
    return res;
}

int last[mo], next[mo];
LL t[mo], v[mo];
int l;

void add(int x, LL y, LL z)
{
    t[++l] = y;
    next[l] = last[x];
    last[x] = l;
    v[l] = z;
}

LL SumPhi(LL x)
{
    if(x < N) return phi[x];
    int k = x % mo;
    for(int i = last[k]; i; i = next[i])
        if(t[i] == x) return v[i];
    LL res = (x % mod) * ((x + 1) % mod) % mod * _2 % mod;
    LL r;
    for(LL i = 2; i <= x; i = r + 1)
    {
        LL tmp = x / i;
        r = x / tmp;
        res = ((res - (SumPhi(tmp) * (r - i + 1) % mod)) % mod + mod) % mod;
    }
    add(k, x, res);
    return res;
}

LL SumPow(LL R, LL L, int k)
{
    LL r = SumOfPow(R, k);
    LL l = SumOfPow(L, k);
    return ((r - l) % mod + mod) % mod;
}

LL solve(LL n, int k)
{
    LL res = 0;
    LL r;
    for(LL i = 1; i <= n; i = r + 1)
    {
        LL tmp = n / i;
        r = n / tmp;
        res = (res + SumPow(r, i - 1, k) * SumPhi(tmp)  % mod) % mod;
    }

    res = res * 2 % mod;

    res = ((res - SumOfPow(n, k)) % mod + mod) % mod;
    return res;
}

int main()
{
  //  freopen("in.txt", "r", stdin);
    init();
    LL n;
    int k;
    scanf("%lld%d", &n, &k);
    printf("%lld\n", solve(n, k));
    return 0;
}