BZOJ2694: Lcm

题目：http://www.lydsy.com/JudgeOnline/problem.php?id=2694

题解：令f[i]表示i是否有平方因子，则f[i]是积性函数，mu[i]表示莫比乌斯函数。

经过balabala的推导，我们得出ans=sigma（f[i/j]*mu[j]*j*j*sum(n/i,m/i))）  sum(x,y)=x*(x+1)/2*y*(y+1)/2

然后我们定义新函数 g[i]=sigma(f[i/d]*mu[d]*d*d) 因为积性函数的狄利克雷卷积仍然是积性函数，所以我们考虑把g数组线筛出来，然后就可以做到sqrt（n）回答询问了。

考虑i%p[j]==0的这部分（初值和i%p[j]!=0可以很简单算出来），如果k=i*p[j]中有p[j]的次数超过2，那么g[k]=0

这是因为我们要在f 和 mu 上分 p[j]的指数，>2时由鸽巢原理知必有一个分到2个以上，那么乘积就是0.

否则 p[j] 的指数为2，我们必须在 f 上分一个，mu上分一个，这样g[k]=g[i/p[j]]*-p[j]*p[j]*p[j]   （第一个p[j]是分到f上的，负号是给 mu 的，p[j]*p[j]则是d*d，还是利用了积性函数的性质）

既然是积性函数并且i/p[j]和p[j]*p[j]互质，那么g[k]就等于g[t]*g[p[j]*p[j]]  注意t==1时要特判。

然后这题就做完了。

因为模数奇特所以直接爆int即可。

代码：

#include<cstdio>

#include<cstdlib>

#include<cmath>

#include<cstring>

#include<algorithm>

#include<iostream>

#include<vector>

#include<map>

#include<set>

#include<queue>

#include<string>

#define inf 1000000000

#define maxn 4000000+5

#define maxm 4000000

#define eps 1e-10

#define ll long long

#define pa pair<int,int>

#define for0(i,n) for(int i=0;i<=(n);i++)

#define for1(i,n) for(int i=1;i<=(n);i++)

#define for2(i,x,y) for(int i=(x);i<=(y);i++)

#define for3(i,x,y) for(int i=(x);i>=(y);i--)

#define for4(i,x) for(int i=head[x],y=e[i].go;i;i=e[i].next,y=e[i].go)

#define for5(n,m) for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)

#define mod 1073741823

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=10*x+ch-'0';ch=getchar();}

    return x*f;

}
int tot,p[maxn],g[maxn];
bool v[maxn];
void get()
{
    g[1]=1;
    for2(i,2,maxm)
    {
        if(!v[i])p[++tot]=i,g[i]=i-i*i;
        for1(j,tot)
        {
            int k=i*p[j];
            if(k>maxm)break;
            v[k]=1;
            if(i%p[j])g[k]=g[i]*g[p[j]];
            else 
            {
                int t=i/p[j];
                if(t%p[j]==0)g[k]=0;
                else g[k]=-g[t]*p[j]*p[j]*p[j];
                break;
            }
        }
    }
    for1(i,maxm)g[i]+=g[i-1];
}
inline int sum(int n,int m)
{
    return n*(n+1)*m*(m+1)/4;
}

int main()

{

    freopen("input.txt","r",stdin);

    freopen("output.txt","w",stdout);
    get();

    int T=read();
    while(T--)
    {
        int n=read(),m=read(),ans=0;
        if(n>m)swap(n,m);
        for(int i=1,j;i<=n;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans+=sum(n/i,m/i)*(g[j]-g[i-1]);
        }
        printf("%d\n",ans&mod);
    }

    return 0;

}