zoukankan      html  css  js  c++  java
  • BZOJ 3601: 一个人的数论

    小清新数学题,全是套路真是舒适无比的说

    我们先化一下式子:

    [ans=sum_{i=1}^n [gcd(i,n)==1] i^d ]

    [=sum_{i=1}^n sum_{t|i,t|n} mu(t) imes i^d ]

    [=sum_{t|n}mu(t) sum_{i=1}^{frac{n}{t}} (it)^d ]

    [=sum_{t|n}mu(t) imes t^d imes sum_{i=1}^{frac{n}{t}} i^d ]

    后面的那个(sum_{i=1}^{frac{n}{t}} i^d)我们都不要太熟悉的说,自然数幂和

    然后我们开始冷静分析一下,(n)很大必然是不能直接求的,也不能直接枚举(n)的约数因为实在太多了

    看起来GG了,我们反过来审视一下题目,发现它给我们的是(n)的质因数分解,看到这个我们再回头看一下式子

    (mu(t))积性函数(t^d)完全积性函数,那么如果后面的自然数幂和换成一个积性函数的话我们就可以对每一个(p_i)单独计算出答案然后乘起来即可

    拉格朗日插值中我们得知(sum_{i=1}^n i^d)是一个关于(n)(d+1)次多项式,如果我们假设这个多项式的系数是(a)的话,则有:

    [sum_{i=1}^n i^d=sum_{i=1}^{d+1} a_i imes n^i ]

    那么假设我们知道了(a),原式子就变成了:

    [=sum_{i=1}^{d+1} a_i imes sum_{t|n}mu(t) imes t^d imes (frac{n}{t})^i ]

    后面的那个是什么,狄利克雷卷积啊!我们知道它也是个积性函数,然后就可以算了

    等下你说就算是(p_i^{alpha_i})的约数也有(alpha_i)个啊,不能一个个枚举。想一下(mu(t))的定义吧,当(t)有平方质因子是(mu(t)=0)

    那么也就意味着我们只用考虑(t=1)(t=p_i)两种情况了,顺利解决

    最后对于那个(a)的计算显然可以用拉格朗日插值的插系数的(O(n^2))写法(不过需要手写多项式除法,比较麻烦),我们看到数据范围(dle 100)直接暴力上高斯消元就可以了

    #include<cstdio>
    #include<iostream>
    #define RI register int
    #define CI const int&
    using namespace std;
    const int N=1005,mod=1e9+7;
    int w,d,p[N],alp[N],a[N][N],val[N],ans;
    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;
    }
    inline int quick_pow(int x,int p=mod-2,int mul=1)
    {
        for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
    }
    inline void Gauss(CI n)
    {
        RI i,j,k; for (i=1;i<=n;++i)
        {
            for (j=i;!a[j][i];++j); swap(a[i],a[j]); swap(val[i],val[j]);
            for (j=i+1;j<=n;++j)
            {
                int tp=1LL*a[j][i]*quick_pow(a[i][i])%mod; dec(val[j],1LL*val[i]*tp%mod);
                for (k=i;k<=n;++k) dec(a[j][k],1LL*a[i][k]*tp%mod);
            }
        }
        for (i=n;i;--i)
        {
            for (j=i+1;j<=n;++j) dec(val[i],1LL*val[j]*a[i][j]%mod);
            val[i]=1LL*val[i]*quick_pow(a[i][i])%mod;
        }
    }
    inline int calc(CI x,int ret=1)
    {
        for (RI i=1;i<=w;++i)
        {
            int tp=quick_pow(p[i],1LL*alp[i]*x%(mod-1));
            if (alp[i]) dec(tp,1LL*quick_pow(p[i],d)*quick_pow(p[i],1LL*(alp[i]-1)*x%(mod-1))%mod);
            ret=1LL*ret*tp%mod;
        }
        return ret;
    }
    int main()
    {
        RI i,j; for (scanf("%d%d",&d,&w),i=1;i<=w;++i) scanf("%d%d",&p[i],&alp[i]);
        for (i=1;i<=d+1;++i) for (j=1;j<=d+1;++j)
        a[i][j]=quick_pow(i,j),j<=i&&(inc(val[i],quick_pow(j,d)),0);
        for (Gauss(d+1),i=1;i<=d+1;++i) inc(ans,1LL*val[i]*calc(i)%mod);
        return printf("%d",ans),0;
    }
    
  • 相关阅读:
    大数运算(涉及到格式问题)
    UltraEdit
    汉化eclipse3.6.2
    安装Microsoft SQL Server Management Studio Express是报错29506
    Java相对路径/绝对路径
    恢复Unbuntu的启动项
    UNC路径
    make: g++:命令未找到
    找到个学习html的网站
    HDU 3756 三分
  • 原文地址:https://www.cnblogs.com/cjjsb/p/12252461.html
Copyright © 2011-2022 走看看