EOJ Monthly 2019.11 E. 数学题

最新推荐文章于 2021-02-07 11:31:59 发布

xufour

最新推荐文章于 2021-02-07 11:31:59 发布

阅读量376

点赞数

分类专栏：反演文章标签：数论

本文链接：https://blog.csdn.net/qq_43206953/article/details/103323444

版权

反演专栏收录该内容

2 篇文章 0 订阅

订阅专栏

链接EOJ Monthly 2019.11 E. 数学题

题意

单点时限: 4.0 sec

内存限制: 256 MB

“飞机已经降落虹桥机场，还将滑行一段时间，请保持安全带的状态。欢迎您再次乘坐东方航空公司的班机，下次旅途再见。”

飞机还在跑道上滑行，到航站楼还有不少距离。Cuber QQ 开始思考他在银川没有做出来的题目。

题目是这样的：统计 k 元组个数 (a1,…,ak),1≤ai≤n 使得 gcd(a1,…,ak,n)=1。

显然当 k=1 时，答案就是欧拉函数。

Cuber QQ 为了方便解答，定义 f(n,k) 为满足要求的 k 元组个数，现在需要求出 ∑i=1nf(i,k)，由于结果可能很大，所以需要对 998 244 353 取模。

Cuber QQ 打算在下飞机前思考出来。但好像还是不会。

输入格式
输入数据只有一行，包含两个函数 n,k(1≤n≤109,1≤k≤1000) 。

输出格式
输出只包含一行，表示求和的结果。

样例
input
4 2
output
24
——————————————————————————————

题意

求 $\sum_{i=1}^{n}\sum_{a_1=1}^{i}\sum_{a_2=1}^{i}.....\sum_{a_k}^{i}[(a_1,a_2...a_k,i)=1]$

解析：

莫比乌斯反演 + 杜教筛 + 伯努利数求幂方和
我们先后面这一段求
$\sum_{i=1}^{n}f$
$f(n)=\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}.....\sum_{a_k}^{n}[(a_1,a_2...a_k,n)=1]$
将等式变为
$\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}.....\sum_{a_k}^{n}[(a_1,a_2...a_k,n)=1] = \sum_{a_1=1}^{n}\sum_{a_2=1}^{n}.....\sum_{a_k}^{n} \sum_{d}^{n}\mu(d) (d|a_1，d|a_2.....d|a_k，d|n)$

枚举d 得
$\sum_{d|n}^{n}\mu(d)\sum_{a_1=1}^{\frac{n}{d}}\sum_{a_2=1}^{\frac{n}{d}}.....\sum_{a_k=1}^{\frac{n}{d}} 1 = \sum_{d|n}^{n}\mu(d)[\frac{n}{d}]^k$
反演一下就是
$\sum_{d|n}^{n}\mu(d)[\frac{n}{d}]^k = \sum_{d|n}^{n}\mu(\frac{n}{d}) n^k = n^k\sum_{d|n}^{n}\mu(\frac{n}{d})$

之后就是杜教筛的环节
原来的函数可以以化简成 $n^k * \mu$
令 $g = I$
$n^k *\mu *I = n^k * e$
杜教筛得
$\sum g*f - \sum_{d=2}^n g(d) *S(\frac{n}{d}) =\sum_{i=1}^{n} i^k - \sum_{d=2}^n g(d) *S(\frac{n}{d})$

之后就是没什么意思的分块整除杜教筛环节了
前面的幂次方可以用伯努利数预处理
整体时间复杂度为 $(n^\frac{2}{3} + k\sqrt n)$

代码

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

#define ll long long
const int MAX = 1e7+9;
const int MAXN= 1007 ;
const ll MOD=998244353;

int pri[MAX+10],tot,k,n;
bool zs[MAX+10];
ll phi[MAX+10];
ll C[1007][1007],inv[1007],B[1007],temp[1007];
ll qpow(ll a,ll b){
    ll e=1;
    while(b){
        if(b%2)e=(e*a)%MOD;
        a=(a*a)%MOD;
        b=b/2;
    }
    return e;
}
void init()//伯努利处理
{
    //预处理组合数
    C[0][0] = 1;
    for (int i = 1; i < MAXN; i++)
    {
        C[i][0] = 1;
        for (int j = 1; j <= i; j++)
        {
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
        }
    }
    //预处理逆元
    inv[1] = 1;
    for (int i = 2; i < MAXN; i++)
    {
        inv[i] = (ll)inv[MOD % i] * (MOD - MOD / i) % MOD;
    }
    //预处理伯努利数
    B[0] = 1;
    for (int i = 1; i < MAXN; i++)
    {
        B[i] = 0;
        for (int k = 0; k < i; k++)
        {
            B[i] = (B[i] + (ll)C[i + 1][k] * B[k] % MOD) % MOD;
        }
        B[i] = ((ll)B[i] * (-inv[i + 1]) % MOD + MOD) % MOD;
    }
}

void pre()//欧拉筛
{
    zs[1]=true;phi[1]=1;
    for(int i=2;i<=MAX;++i)
    {
        if(!zs[i])
            pri[++tot]=i,phi[i]=qpow(i,k)-1;
        for(int j=1;j<=tot&&pri[j]*i<=MAX;++j)
        {
            zs[i*pri[j]]=true;
            if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]]%MOD;
            else{phi[i*pri[j]]=phi[i]*qpow(pri[j],k)%MOD;break;}
        }
    }
    for(int i=1;i<=MAX;++i)phi[i]=(phi[i] + phi[i-1])%MOD;

}
ll SSS(ll x)// 求幂方和
{
        x++;
        x %= MOD;
        ll tmp = x, ans = 0;
        for (int i = 1; i <= k + 1; i++)
        {
            ans = (ans + (ll)C[k + 1][i] * B[k + 1 - i] % MOD * x % MOD) % MOD;
            x = (ll)x * tmp % MOD;
        }
        ans = (ll)ans * inv[k + 1] % MOD;
        return ans ;

}
map<ll,ll> M;
ll Solve(ll x)//杜教筛
{
    if(x<=MAX)return phi[x];
    if(M[x])return M[x];
    ll ret=0;
    for(ll i=2,j;i<=x;i=j+1)
    {
        j=x/(x/i);
        ret=(ret+(j-i+1)*Solve(x/i))%MOD;
    }

    return M[x]=(SSS(x)-ret+MOD)%MOD;
}



int main()
{

    cin>>n>>k;
    init();
    pre();
    printf("%lld\n",Solve(n));
    return 0;
}