zoukankan      html  css  js  c++  java
  • bzoj3512: DZY Loves Math IV

    开始莫反强行拿掉欧拉函数搞出了一个这样的柿子。。。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
    ",ans);
        
        return 0;
    }
  • 相关阅读:
    OAuth认证介绍及腾讯微博OAuth认证示例
    Android PopupWindow介绍及实现菜单效果
    Android 输入法键盘和activity页面遮挡问题解决
    eclipse中android项目的编译过程分析
    Android Tween动画之RotateAnimation实现图片不停旋转
    eclipse 文件同步插件
    关于移动Web应用程序开发 HTML5、高性能JavaScript篇、Css的几篇较好博客
    Android 记录和恢复ListView滚动的位置的三种方法
    Apache的功能模块
    如何防止自己网站的图片被其他网站所盗用,从而导致自己网站流量的损失【apache篇】
  • 原文地址:https://www.cnblogs.com/AKCqhzdy/p/10592752.html
Copyright © 2011-2022 走看看