传送门: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 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;
}