BZOJ3512：DZY Loves Math IV（杜教筛）

最新推荐文章于 2020-09-23 16:08:09 发布

Fat_tu

最新推荐文章于 2020-09-23 16:08:09 发布

阅读量640

点赞数

CC 4.0 BY-SA版权

分类专栏：数论杜教筛文章标签：数论

本文链接：https://blog.youkuaiyun.com/q582116859/article/details/78931750

数论同时被 2 个专栏收录

19 篇文章

订阅专栏

杜教筛

3 篇文章

订阅专栏

一道涉及数论的题目，使用杜教筛方法解决，复杂度为nm的平方根，通过巧妙的记忆化实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

很巧妙的一道题
题面
题意：

\sum i = 1 n \sum j = 1 m φ (i j)

$\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)$
n≤1e5，m≤1e9

显然杜教筛是这个问题的子问题，而我又没学过除了杜教筛以外的算法
所以这题就是杜教筛了。

我在语文课上突然想起了markdown写公式的恐怖。
这里写图片描述
这个试子复杂度上界为 $n\sqrt{m}$ ，玄学记忆化就可以过。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <map>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=1001000,nn=1000000,mo=1e9+7;

int n,m,ans;
int prime[N],num,phi[N],w[N];
LL sphi[N];
bool b[N];
map<int ,LL> mp;
map<LL,LL> mmp;

LL getp(int x)
{
    if(x<=nn)
    return sphi[x];
    if(mp.count(x))
    return mp[x];
    LL res=(LL)(x+1)*x/2;
    for(LL i=2,last;i<=x;i=last+1)
    {
        last=x/(x/i);
        res-=(last-i+1)*getp(x/i);
    }
    mp[x]=res;
    return res;
}

int work(int x,int y)
{
    if(mmp.count(1000000000ll*x+y))
    return mmp[1000000000ll*x+y];

    if(x==1)
    return getp(y)%mo;
    if(y==0)
    return 0;
    if(y==1)
    return phi[x];
    int hy=x/w[x],tu=sqrt(w[x]),res=0;
    for(int i=1;i<=tu;i++)
    if(w[x]%i==0)
    {
        int oi=w[x]/i;
        res=(res+(LL)phi[oi]*work(i,y/i)%mo)%mo;
        if(i*i!=w[x])
        res=(res+(LL)phi[i]*work(oi,y/oi)%mo)%mo;
    }
    res=(LL)res*hy%mo;
    mmp[1000000000ll*x+y]=res;
    return res;
}

int main()
{
    phi[1]=1;
    for(int i=2;i<=nn;i++)
    {
        if(!b[i])
        {
            prime[++num]=i;
            phi[i]=i-1;
            w[i]=i;
        }
        for(int j=1;j<=num&&prime[j]*i<=nn;j++)
        {
            int k=i*prime[j];
            b[k]=1;
            if(i%prime[j]==0)
            {
                phi[k]=phi[i]*prime[j];
                w[k]=w[i];
                break;
            }
            phi[k]=phi[i]*(prime[j]-1);
            w[k]=w[i]*prime[j];
        }
    }

    for(int i=1;i<=nn;i++)
    sphi[i]=sphi[i-1]+phi[i];

    cin>>n>>m;
    for(int i=n;i>=1;i--)
    ans=(ans+work(i,m))%mo;

    cout<<ans<<endl;

    return 0;
}