• BZOJ 3512: DZY Loves Math IV


    这是一个悲伤的故事,我之前递归求(S(n,m))的时候忘记给记忆化数组赋值了,然后就跑得很慢(废话)

    然后我一直以为自己的杜教筛写的太辣鸡了,分解质因数太辣鸡了,白调了2h……(自闭ing)

    总的来说这题确实是妙,又教会我一个常用(?)套路的说。

    首先我们注意到(n)范围不大,因此我们考虑枚举(n),然后设一个(S(n,m)=sum_{i=1}^m phi(ni)),考虑到(phi)的性质,我们可以把(n)的一些重复的质因数拿出来,形式化的:

    (n=prod_{i} p_i^{a_i}),那么(phi(n)=phi(prod_i p_i) imes prod_i p_i^{a_i-1}),那么我们令(w=prod_i p_i,y=prod_i p_i^{a_i-1}),则:

    [S(n,m)=y imes sum_{i=1}^m phi(wi)\=y imes sum_{i=1}^m phi(frac{w}{gcd(w,i)}) imes phi(i) imes gcd(i,w)\=y imes sum_{i=1}^m phi(frac{w}{gcd(w,i)}) imes phi(i) imes sum_{d|gcd(w,i)}phi(d)\=y imes sum_{i=1}^m phi(i) imes sum_{d|gcd(w,i)}phi(frac{wd}{gcd(w,i)})\=y imes sum_{i=1}^m phi(i) imes sum_{d|gcd(w,i)}phi(frac{w}{d})\=y imes sum_{i=1}^m phi(i) imes sum_{d|w,d|i}phi(frac{w}{d})\=y imes sum_{d|w} phi(frac{w}{d}) imes sum_{i=1}^{lfloorfrac{m}{d} floor}phi(di)\=y imes sum_{d|w} phi(frac{w}{d}) imes S(d,lfloorfrac{m}{d} floor) ]

    中间的代换有一步很妙,运用(n=sum_{d|n} phi(d))进行变换,把(gcd)给去掉了,看来这个技巧和(mu)(gcd)都要掌握啊

    然后我们得到了一个递归式,当(n=1)(S(n,m)=sum_{i=1}^m phi(i)),用杜教筛筛出来即可,同时我们将(S(n,m))也记忆化,来算一下复杂度

    首先杜教筛的复杂度是(O(m^{frac{2}{3}})),并且是针对全局的,然后每一个(S(n,m))枚举分解(n)的复杂度是(O(sqrt n))的,然后每一个(n)的对应的(m)的取值只用(sqrt m)种,因此总复杂度是(O(n(sqrt n+sqrt m)+m^{frac{2}{3}}))

    PS:以下的代码用了优化的不用map的杜教筛,比起map还是快了不少

    #include<cstdio>
    #include<map>
    #include<cmath>
    #define RI register int
    #define CI const int&
    using namespace std;
    const int N=100005,M=1e6,mod=1e9+7;
    int n,m,sum_phi[M+5],Phi[M+5],phi[M+5],prime[M+5],cnt,ans,d,id1[M+5],id2[M+5],idx;
    bool vis[M+5]; map <int,int> s[N];
    inline void inc(int &x,CI y)
    {
        if ((x+=y)>=mod) x-=mod;
    }
    inline void dec(int& x,CI y)
    {
        if ((x-=y)<0) x+=mod;
    }
    #define P(x) prime[x]
    inline void init(CI n)
    {
        RI i,j; for (phi[1]=vis[1]=1,i=2;i<=n;++i)
        {
            if (!vis[i]) prime[++cnt]=i,phi[i]=i-1;
            for (j=1;j<=cnt&&i*P(j)<=n;++j)
            {
                vis[i*P(j)]=1; if (i%P(j)) phi[i*P(j)]=1LL*phi[i]*(P(j)-1)%mod;
                else phi[i*P(j)]=1LL*phi[i]*P(j)%mod;
            }
        }
        for (i=1;i<=n;++i) sum_phi[i]=sum_phi[i-1],inc(sum_phi[i],phi[i]);
    }
    #define ID(x) (x<=d?id1[x]:id2[m/x])
    inline int get_phi(CI x)
    {
        if (x<=M) return sum_phi[x]; if (Phi[ID(x)]) return Phi[ID(x)];
        int ret=(1LL*x*(x+1)/2LL)%mod; for (RI l=2,r;l<=x;l=r+1)
        r=x/(x/l),dec(ret,1LL*(r-l+1)*get_phi(x/l)%mod); return Phi[ID(x)]=ret;
    }
    #undef ID
    inline int S(int n,CI m)
    {
        if (n==1) return get_phi(m); if (m==1) return phi[n]; if (!m) return 0;
        if (s[n].count(m)) return s[n][m]; RI i; int w=1,y=1,ret=0,t=n;
        for (i=1;i<=cnt&&P(i)*P(i)<=n;++i) if (n%P(i)==0)
        for (w*=P(i),n/=P(i);n%P(i)==0;y*=P(i),n/=P(i));
        if (n>1) w*=n; for (i=1;i*i<=w;++i) if (w%i==0)
        {
            inc(ret,1LL*phi[w/i]*S(i,m/i)%mod);
            if (i*i!=w) inc(ret,1LL*phi[i]*S(w/i,m/(w/i))%mod);
        }
        return s[t][m]=1LL*y*ret%mod;
    }
    #undef P 
    int main()
    {
        scanf("%d%d",&n,&m); init(M); d=sqrt(m); for (RI l=1,r;l<=m;l=r+1)
        r=m/(m/l),(m/l<=d)?id1[m/l]=++idx:id2[r]=++idx;
        for (RI i=1;i<=n;++i) inc(ans,S(i,m)); return printf("%d",ans),0;
    }
    
  • 相关阅读:
    雅虎天气API调用
    HttpOperater
    HttpOperater-模拟HTTP操作类
    页面局部加载,适合Ajax Loading场景(Demo整理)
    FTPHelper-封装FTP的相关操作
    使用SOCKET实现TCP/IP协议的通讯
    IIS目录禁止执行权限
    Oracle10g 安装步骤
    SQL Server 2008、SQL Server 2008R2 自动备份数据库
    SQL列转行
  • 原文地址:https://www.cnblogs.com/cjjsb/p/12257795.html
Copyright © 2020-2023  润新知