OI回忆录

【数论】bzoj4174tty的求助

题目链接
(ps:这是一道权限题…)

计算

n=1Nm=1Mk=0m1nk+xm

其中NM<=5000000<x<=100000x精确到小数点后8
位小数。

根本不会的说…
然后去看po姐博客就一脸懵逼辣…

以下摘自popoqqq博客,并加上蒟蒻的辣鸡理解…
单独考虑最后一个和式

k=0m1nk+xm

=k=0m1(nk(mod)m+xm+nkmnk(mod)mm)

上式是考虑nkm的倍数部分。

然后分开考虑每一个独立的和式。
d=gcd(n,m),这个变量将贯穿整篇博客。

k=0m1nk(mod)m+xm

这相当于在枚举m的剩余系,然后很明显剩余系重复了d次,所以只需要枚举一个剩余系即可。
=dk=0md1kd+xm

和之前同样的考虑方式,考虑xm的倍数部分。
=d(mdxx(mod)mm+k=0md1kd+x(mod)mm)

很明显,kd+x(mod)m是不会比2m大的…
=d(xx(mod)md+k=0md1[kd+x(mod)mm])

移项
=d(xx(mod)md+k=0md1[kmx(mod)md])

=d(xx(mod)md+mdmx(mod)md)

=d(xx(mod)md+x(mod)md)

=dxd

对于第二个和式:

k=0m1nkm

根据和式的性质,可以将常数提取出来。
=nk=0m1km=nm(m1)2m

=nmn2

对于第三个和式,化简的过程与第一个和式的部分相似。

k=0m1nk(mod)mm

这相当于在枚举m的剩余系,然后很明显剩余系重复了d次,所以只需要枚举一个剩余系即可。
=dk=0md1kdm=d2mmd(md1)2

=md2

接着,我们来整理一下这些柿子:

n=1Nm=1M(dxd+nmn2md2)

=12n=1Nm=1M(2dxd+nmnm+d)

把和式展开,其中Sum(i)表示从1i的公差为1的等差数列的前n项和,即Sum(i)=i(i+1)2
=12(Sum(N)Sum(M)NSum(M)MSum(N)+n=1Nm=1M(2dxd+d))

对于最后的和式中的d=gcd(n,m),我们考虑每一个d出现的次数,这样就可以计算每一个相同的d对答案的贡献。
=12(Sum(N)Sum(M)NSum(M)MSum(N)+

d=1min(N,M)(2dxd+d)g=1min(nd,md)μ(g)NdgMdg))

然后我们就枚举d,然后用经典的莫比乌斯反演求一下gcd(i,j)==d(1<=i<=N,1<=j<=M)的个数

#include <iostream>
#include <cstdio>
#define LL long long int
#define mod 998244353
#define MAXN 500005
using namespace std;

int n, m, prime[MAXN], tot;
LL u[MAXN];
bool flag[MAXN];
double x;

inline void init()
{
    u[1]=1;
    for(int i=2, j, k;i<MAXN;++i)
    {
        if(!flag[i])prime[++tot]=i, u[i]=-1;
        for(j=1;j<=tot&&(k=i*prime[j])<MAXN;++j)
        {
            flag[k]=1;
            if(i%prime[j]==0){u[k]=0;break;}
            u[k]=-u[i];
        }
    }
    for(int i=1;i<MAXN;++i)
        u[i]=(u[i]+u[i-1]+mod+mod)%mod;
}

LL sum(LL n){return n*(n+1)/2%mod;}

int main()
{
    init();
    scanf("%d%d%lf",&n,&m,&x);
    LL ans=((sum(n)*sum(m)-n*sum(m)-m*sum(n))%mod+mod)%mod, k;
    if(n>m)swap(n,m);
    for(int i=1, nn, mm, last;i<=n;++i)
    {
        k=i+int(x/i)*i*2;
        nn=n/i, mm=m/i;
        for(int j=1;j<=nn;j=last+1)
        {
            last=min(nn/(nn/j),mm/(mm/j));
            ans=(ans+k*(u[last]-u[j-1])%mod*(nn/j)%mod*(mm/j))%mod;
        }
    }
    cout<<ans*499122177%mod<<'\n';
    return 0;
}
阅读更多
版权声明:转载请注明出处,O(∩_∩)O谢谢 https://blog.csdn.net/acmer_cheer/article/details/51283073
文章标签: 数论 bzoj4174
个人分类: 数论
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭