BZOJ3561 - DZY Loves Math VI

最新推荐文章于 2021-09-28 08:45:10 发布

转载最新推荐文章于 2021-09-28 08:45:10 发布 · 98 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/VisJiao/p/BZOJ3561.html

本文探讨了一道数学竞赛题目，目标是求解给定范围内所有数对的最小公倍数与最大公约数幂次的总和。通过巧妙的数学转换和算法优化，将复杂度降低到可接受范围，适用于大规模数据集。代码实现采用C++，利用预处理和枚举策略高效解决该问题。

Portal

Description

共\(T(T\leq3)\)组测试数据。给出\(n,m(n\leq5\times10^5)\)，求
\[ \sum_{i=1}^n \sum_{j=1}^m lcm(i,j)^{gcd(i,j)}\]

Solution

喜闻乐见推推推。
\[\begin{align*} ans &= \sum_{i=1}^n \sum_{j=1}^m lcm(i,j)^{gcd(i,j)} \\ &= \sum_{d=1}^{+∞} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d](\frac{ij}{d})^d \\ &= \sum_{d=1}^{+∞} \sum_{i=1}^{⌊\frac{n}{d}⌋} \sum_{j=1}^{⌊\frac{m}{d}⌋} [gcd(i,j)=1](ijd)^d \\ &= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1) \sum_{d_1|i}^{⌊\frac{n}{d}⌋} \sum_{d_1|j}^{⌊\frac{m}{d}⌋} (ij)^d \\ &= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1) \sum_{i=1}^{⌊\frac{n}{dd_1}⌋} \sum_{j=1}^{⌊\frac{m}{dd_1}⌋} (id_1\cdot jd_1)^d \\ &= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1)d_1^{2d} \sum_{i=1}^{⌊\frac{n}{dd_1}⌋} i^d\sum_{j=1}^{⌊\frac{m}{dd_1}⌋} j^d \\ \end{align*}\] 于是我们枚举\(d\)，预处理出\(f(x)=\sum_{i=1}^x i^d\)，然后枚举\(d_1\)。

根据调和级数总复杂度为\(O(nlogn)\)。

Code

//DZY Loves Math VI
#include <algorithm>
#include <cstdio>
using std::swap;
typedef long long lint;
const int N=5e5+10;
const int P=1e9+7;
int prCnt,pr[N]; bool prNot[N];
int mu[N];
void init(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!prNot[i]) pr[++prCnt]=i,mu[i]=-1;
        for(int j=1;j<=prCnt;j++)
        {
            int x=i*pr[j]; if(x>n) break;
            prNot[x]=true;
            if(i%pr[j]) mu[x]=-mu[i]; else break;
        }
    }
}
int powD[N],S[N];
int pow(int x,int y)
{
    lint r=1,t=x;
    for(int i=y;i;i>>=1,t=(t*t)%P) if(i&1) r=(r*t)%P;
    return r;
}
int main()
{
    int task=1; init(5e5);
    while(task--)
    {
        int n,m; scanf("%d%d",&n,&m);
        if(n>m) swap(n,m);
        int ans=0;
        for(int i=1;i<=m;i++) powD[i]=1;
        for(int d=1;d<=n;d++)
        {
            int k1=n/d,k2=m/d,ansD=0;
            for(int i=1;i<=k2;i++) powD[i]=(1LL*powD[i]*i)%P;
            for(int i=1;i<=k2;i++) S[i]=(S[i-1]+powD[i])%P;
            for(int d1=1;d1<=k1;d1++) 
            {
                lint t=1LL*(mu[d1]+P)%P*powD[d1]%P*powD[d1]%P;
                t=t*S[k1/d1]%P*S[k2/d1]%P;
                ansD=(ansD+t)%P;
            }
            ans=(ans+1LL*ansD*pow(d,d))%P;
        }
        printf("%d\n",ans);
    }
    return 0;
}

转载于:https://www.cnblogs.com/VisJiao/p/BZOJ3561.html