DZY Loves Math VI

Posted 2020-12-30 wolfycz

Description
给定正整数n,m。求
[sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j)^{gcd(i,j)}]

Input
一行两个整数n,m。

Output
一个整数，为答案模1000000007后的值。

Sample Input
5 4

Sample Output
424

HINT
数据规模：
1<=n,m<=500000，共有3组数据。

首先推柿子
[sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j)^{gcd(i,j)}]
[sumlimits_{i=1}^nsumlimits_{j=1}^m(dfrac{i imes j}{gcd(i,j)})^{gcd(i,j)}]
[sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}(dfrac{d^2 imes i imes j}{d})^d[gcd(i,j)==1]]
[sumlimits_{d=1}^nd^dsumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}i^d imes j^dsumlimits_{x|i,x|j}mu(x)]
[sumlimits_{d=1}^nd^dsumlimits_{x=1}^{lfloorfrac{n}{d} floor}mu(x) imes x^{2d}sumlimits_{i=1}^{lfloorfrac{n}{dx} floor}i^dsumlimits_{j=1}^{lfloorfrac{m}{dx} floor}j^d]
然后设(T=dx)，然后就。。。布星，推不下去了
发现m<=5e5，那不就得了，直接这样求就好，枚举d的时候，维护一下(f(n)=sumlimits_{i=1}^n i^d)即可，复杂度约为(O(mlog m))

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
    int x=0,f=1;char ch=getchar();
    for (;ch<‘0‘||ch>‘9‘;ch=getchar())  if (ch==‘-‘)    f=-1;
    for (;ch>=‘0‘&&ch<=‘9‘;ch=getchar())    x=(x<<1)+(x<<3)+ch-‘0‘;
    return x*f;
}
inline void print(int x){
    if (x>=10)  print(x/10);
    putchar(x%10+‘0‘);
}
const int N=5e5,p=1e9+7;
int prime[N+10],miu[N+10];
bool inprime[N+10];
void prepare(){
    int tot=0; miu[1]=1;
    for (int i=2;i<=N;i++){
        if (!inprime[i])    prime[++tot]=i,miu[i]=-1;
        for (int j=1;j<=tot&&prime[j]*i<=N;j++){
            inprime[i*prime[j]]=1;
            if (i%prime[j]==0)  break;
            miu[i*prime[j]]=-miu[i];
        }
    }
}
int mlt(int a,int b){
    int res=1;
    for (;b;b>>=1,a=1ll*a*a%p)  if (b&1)    res=1ll*res*a%p;
    return res;
}
int val[N+10],sum[N+10];
int main(){
    prepare();
    int n=read(),m=read(),Ans=0;
    if (n>m)    swap(n,m);
    for (int i=1;i<=m;i++)  val[i]=1;
    for (int d=1;d<=n;d++){
        int x=mlt(d,d),res=0,limn=n/d,limm=m/d;
        for (int i=1;i<=limm;i++)   val[i]=1ll*val[i]*i%p,sum[i]=(sum[i-1]+val[i])%p;
        for (int i=1;i<=limn;i++){
            if (!miu[i])    continue;
            res=(res+1ll*val[i]*val[i]%p*sum[limn/i]%p*sum[limm/i]%p*miu[i]+p)%p;
        }
        Ans=(Ans+1ll*x*res%p)%p;
    }
    printf("%d
",Ans);
    return 0;
}