uestc 811 GCD 杜教筛 + 自然幂和

传送门:UESTC 811

---

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

这题好像不好找题解

题意

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

求∑Ni=1∑Nj=1gcd(i,j)k

---

题解

枚举gcd, 推下式子

ans=∑i=1N∑j=1Ngcd(i,j)k=2∑i=1N∑j=1igcd(i,j)k−∑i=1Ngcd(i,i)k=2∑i=1N∑d|idkϕ(id)−∑i=1Nik=2∑i=1Nik∑d=1⌊Ni⌋ϕ(d)−∑i=1Nik

然后就可以分块做了

---

快速求自然幂和

∑Ni=1ik=1k+1∑k+1i=1Cik+1Bk+1−i(n+1)k
预处理出伯努利数和组合数就可以在O(k+1)时间内求出自然幂和,
k最大为5

---

伯努利数求法:
B0=1,Bk=−1k+1∑k−1i=0Cik+1Bi
求法源于:
∑ki=0Cik+1Bi=0

---

欧拉函数的杜教筛之前题解里涉及到, 就不多讲了
总复杂度O(n23)

---

code:

#include <bits/stdc++.h>using namespace std;typedef long long LL;#define next Nextconst 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;}