首页 > 其他 > 详细

于神之怒加强版 解题报告

时间:2019-03-17 20:43:48      阅读:166      评论:0      收藏:0      [点我收藏+]

于神之怒加强版

题目描述

给定\(n,m,k\),计算
\[ \sum_{i=1}^n \sum_{j=1}^m \mathrm{gcd}(i,j)^k \]
\(1000000007\)取模的结果

说明

\(T\le 2000,1\le N,M,K\le 5000000\)


注意\(k\)不是每次都给的...

可以推出式子
\[ \sum_{T=1}^{\min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{ab=T}a^k\mu(b) \]
然后预处理一下右边,发现右边是\(f=Id^k*\mu\),所以是积性的

因为\(Id^k\)是完全积性函数,所以可以\(O(n)\)预处理

然后又有\(f(p^c)=(p^c)^k-(p^{c-1})^k\)\(p\)是质数,就可以筛了


Code:

#include <cstdio>
#include <algorithm>
using std::min;
const int N=5e6+1;
const int mod=1000000007;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
#define mul(a,b) (1ll*(a)*(b)%mod)
int qp(int d,int k){int f=1;while(k){if(k&1)f=mul(f,d);d=mul(d,d),k>>=1;}return f;}
int Id[N],b[N],f[N],pri[N],ispri[N],k,cnt;
void init()
{
    Id[1]=b[1]=f[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!ispri[i])
        {
            Id[i]=qp(i,k);
            f[i]=Id[i]-1;
            b[i]=i;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<N;j++)
        {
            int x=pri[j]*i;
            ispri[x]=1;
            Id[x]=mul(Id[i],Id[pri[j]]);
            if(i%pri[j])
            {
                b[x]=pri[j];
                f[x]=mul(f[i],f[pri[j]]);
            }
            else
            {
                b[x]=mul(b[i],pri[j]);
                if(b[i]==i) f[x]=add(Id[x],mod-Id[i]);
                else f[x]=mul(f[i/b[i]],f[b[x]]);
                break;
            }
        }
    }
    for(int i=2;i<N;i++) f[i]=add(f[i],f[i-1]);
}
int main()
{
    int T,n,m;scanf("%d%d",&T,&k);
    init();
    while(T--)
    {
        scanf("%d%d",&n,&m);
        int ans=0;
        if(n>m) std::swap(n,m);
        for(int l=1,r;l<=n;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans=add(ans,mul(n/l,mul(m/l,add(f[r],mod-f[l-1]))));
        }
        printf("%d\n",ans);
    }
    return 0;
}

2019.3.17

于神之怒加强版 解题报告

原文:https://www.cnblogs.com/butterflydew/p/10548464.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!