bzoj3512: DZY Loves Math IV

最新推荐文章于 2023-01-12 16:04:40 发布

转载最新推荐文章于 2023-01-12 16:04:40 发布 · 67 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/AKCqhzdy/p/10592752.html

本文详细解析了一道关于杜教筛算法的题目，通过具体的数学推导和代码实现，介绍了如何利用杜教筛求解特定的数论问题，并提供了一个完整的C++实现示例。

开始莫反强行拿掉欧拉函数搞出了一个这样的柿子。。。sigema(1~n)c sigema(1~n*m/c)d u(d)* (1+n/c)*(n/c)/2 * (1+m/(d/c))*(m/(d/c))/2

自以为很对结果后来发现居然没办法把u给处理出来。。。囧

结果发现是个杜教筛裸题的哥哥题

然后这个n这么小肯定就是给你枚举的，令S(n,m)=sigema(1~m)i phi(n*i)，答案就是sigema(1~n)i S(i,m)

考虑S(n,m)怎么算

对于欧拉函数的值，根据柿子可以理解为每个质因子被加入ci次，第一次贡献为pi-1,之后就是p

把n给拆了，每个质因子只保留1个记作w，原式=(n/w)*sigema(1~m)i phi(w*i)

根据欧拉函数的性质可得sigema(1~m)i phi(w*i)=sigema(1~m) phi(w/gcd(w,i))*phi(i)*gcd(w,i)

把gcd拿掉变成sigema(1~m) phi(i)*phi(w/gcd(w,i))*sigema(d|gcd(w,i))phi(d)

d和w/gcd(w,i)明显互质，压进去就变成 sigema(1~m) phi(i)*sigema(d|gcd(w,i))phi(w/d)

改变枚举顺序，再把整除去掉变成sigema(d|w)phi(w/d)*sigema(1~m/d)i phi(i*d)=sigema(d|w)phi(w/d)*S(d,m/d)

那么n=1直接杜教筛，m<=1特判一下，S(n,m)=(n/w)*sigema(d|w)phi(w/d)*S(d,m/d)，记忆化搜索就过了

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<map>
using namespace std;
typedef long long LL;
const int _=1e2;
const int maxn=1e6+_;
const LL mod=1e9+7;

int pr,prime[maxn];bool v[maxn];
int phi[maxn];LL sphi[maxn];
void getprime()
{
    phi[1]=sphi[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(v[i]==false)
            prime[++pr]=i,phi[i]=i-1;
        for(int j=1;j<=pr&&i*prime[j]<maxn;j++)
        {
            int k=i*prime[j];
            v[k]=true;
            if(i%prime[j]==0)
                {phi[k]=phi[i]*prime[j];break;}
            else phi[k]=phi[i]*(prime[j]-1);
        }
        sphi[i]=sphi[i-1]+phi[i];if(sphi[i]>=mod)sphi[i]-=mod;
    }
}
map<int,int>PHI;
int getphi(int n)
{
    if(n<maxn)return phi[n];
    if(PHI.find(n)!=PHI.end())return PHI[n];
    int ret=n;
    for(int i=2;i*i<=n;i++)
        if(n%i==0)
        {
            ret=ret/i*(i-1);
            while(n%i==0)n/=i;
        }
    if(n!=1)ret=ret/n*(n-1);
    PHI[n]=ret;
    return ret;
}
map<int,LL>SPHI;
LL getsphi(int n)
{
    if(n<maxn)return sphi[n];
    if(SPHI.find(n)!=SPHI.end())return SPHI[n];
    LL ret=(LL)(1+n)*n/2%mod; int last;
    for(int d=2;d<=n;d=last+1)
    {
        last=n/(n/d);
        ret=ret-(last-d+1)*getsphi(n/d)%mod;
        if(ret<0)ret+=mod;
    }
    SPHI[n]=ret;
    return ret;
}

map<pair<int,int>,LL>S;
LL gets(int n,int m)
{
    if(m==0)return 0;
    if(m==1)return getphi(n);
    if(n==1)return getsphi(m);
    pair<int,int> P=make_pair(n,m);
    if(S.find(P)!=S.end())return S[P];
    
    int plen=0,p[10];
    int k=n,w=1;
    for(int i=2;i*i<=k;i++)
        if(k%i==0)
        {
            p[++plen]=i; w*=i;
            while(k%i==0)k/=i;
        }
    if(k!=1)p[++plen]=k, w*=k;
    //debug
    
    LL ret=0;
    int li=(1<<plen)-1;
    for(int zt=li;zt;zt=(zt-1)&li)
    {
        int d=1;
        for(int i=1;i<=plen;i++)
            if(zt&(1<<i-1))d*=p[i];
        ret=ret+getphi(w/d)*gets(d,m/d)%mod;
        if(ret>mod)ret-=mod;
    }
    ret=ret+getphi(w)*gets(1,m)%mod;
    if(ret>mod)ret-=mod;
    ret=(ret*(n/w))%mod;
    
    S[P]=ret;
    return ret;
}

int main()
{
    getprime();
    int n,m;
    scanf("%d%d",&n,&m);
    LL ans=0;
    for(int i=1;i<=n;i++)
    {
        ans=ans+gets(i,m);
        if(ans>=mod)ans-=mod;
    }
    printf("%lld\n",ans);
    
    return 0;
}

转载于:https://www.cnblogs.com/AKCqhzdy/p/10592752.html