[BZOJ4407]于神之怒加强版

BZOJ挂了。。。

先把程序放上来，如果A了在写题解吧。

#include<cstdio>
#include<algorithm>
#define N 5000010
#define ll long long
#define mod (int)(1e9+7)
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
inline int ksm(ll x,ll y)
{
    ll p=1;
    while(y)
    {
        if(y&1)p=p*x%mod;
        x=x*x%mod;y>>=1;
    }
    return p;
}
int prime[N],tot,f[N],p[N],k;
bool vis[N];
void pre(int t)
{
    f[1]=1;
    for(int i=2;i<=t;i++)
    {
        if(!vis[i])
        {
            prime[++tot]=i;
            p[tot]=ksm(i,k);
            f[i]=p[tot]-1;
        }
        for(int j=1;j<=tot&&i*prime[j]<=t;j++)
        {
            vis[i*prime[j]]=1;
            if(!(i%prime[j]))
            {
                f[i*prime[j]]=(ll)f[i]*p[j]%mod;
                break;
            }
            f[i*prime[j]]=(ll)f[i]*f[prime[j]]%mod;
        }
    }
    for(int i=1;i<=t;i++)
    f[i]=(f[i]+f[i-1])%mod;
}
int n[2005],m[2005],t,maxn;
int main()
{
    t=read();k=read();
    for(int i=1;i<=t;i++)
    {
        n[i]=read();
        m[i]=read();
        if(n[i]>m[i])swap(n[i],m[i]);
        maxn=max(maxn,n[i]);
    }
    pre(maxn);
    for(int z=1;z<=t;z++)
    {
        int ans=0;
        for(int i=1,j;i<=n[z];i=j+1)
        {
            j=min(n[z]/(n[z]/i),m[z]/(m[z]/i));
            ans=(ans+((ll)(f[j]+mod-f[i-1])*(n[z]/i)%mod*(m[z]/i)%mod))%mod;
        }
        printf("%d
",ans);
    }
}

update:

$sumlimits_{i=1}^nsumlimits_{j=1}^ngcd(i,j)^k$
$=sumlimits_{d=1}^nd^k*sumlimits_{i=1}^{[frac nd]}sumlimits_{j=1}^{[frac md]}gcd(i,j)==1$
$=sumlimits_{d=1}^nd^k*sumlimits_{i=1}^{[frac nd]}sumlimits_{j=1}^{[frac md]}sumlimits_{t|i,t|j}mu(t)$
$=sumlimits_{d=1}^nd^k*sumlimits_{t=1}^nmu(t)*[frac n{td}]*[frac m{td}]$
令$x=d*t$
原式=$sumlimits_{x=1}^n[frac nx][frac mx]sumlimits_{d|x}mu(frac xd)*d^k$

后面那个显然是积性函数，可以线性筛。线性筛预处理+分块 时间复杂度$O(n+T*sqrt n)$