BZOJ2694：Lcm——包看得懂/看不懂题解

http://www.lydsy.com/JudgeOnline/problem.php?id=2694
Description

对于任意的>1的n gcd(a, b)不是n^2的倍数
也就是说gcd(a, b)没有一个因子的次数>=2
Input
一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M
Output
T行 每行一个整数 表示第i组数据的结果
Sample Input
4
2 4
3 3
6 5
8 3
Sample Output
24
28
233
178
——————————————————————————————————
觉得这一道题很有纪念意义所以准备用markdown好好用公式表达一下。

而且发现我所看的博客都或多或少存在一些问题，所以不再贴他们的博客地址了。

定义：一个函数(f(x))，当x包含平方因子的时候，(f(x)=0)，否则(f(x)=1)

定义：一个函数(sum(x)=sum_{i=1}^x i)

(如果你觉得下面的过程太冗杂了，请跳过大括号区域内的过程)

——————————————————————————————————

首先假设(n<m)，那么我们的答案最初为:

(sum_{i=1}^nsum_{j=1}^m lcm(i,j)*f(gcd(i,j)))

(=sum_{p=1}^nsum_{i=1}^nsum_{j=1}^mfrac{ij}{p}*f(p)*[gcd(i,j)=p])

(=sum_{p=1}^nsum_{i=1}^{lfloorfrac{n}{p} floor}sum_{j=1}^{lfloorfrac{m}{p} floor} pij*f(p)*[gcd(i,j)=1])

(=sum_{p=1}^nsum_{i=1}^{lfloorfrac{n}{p} floor}sum_{j=1}^{lfloorfrac{m}{p} floor} pij*f(p) sum_{d|gcd(i,j)}mu(d))

({)

(=sum_{p=1}^nsum_{i=1}^{lfloorfrac{n}{p} floor}sum_{j=1}^{lfloorfrac{m}{p} floor} pij*f(p) sum_{d|gcd(i,j)}mu(d))

(=sum_{p=1}^nsum_{i=1}^{lfloorfrac{n}{p} floor}sum_{j=1}^{lfloorfrac{m}{p} floor} pij*f(p) sum_{d|iigwedge d|j}mu(d))

(=sum_{p=1}^n p*f(p)sum_{d=1}^{lfloorfrac{n}{p} floor}mu(d)sum_{i=1igwedge d|i}^{lfloorfrac{n}{p} floor}sum_{j=1igwedge d|j}^{lfloorfrac{m}{p} floor} ij)

(=sum_{p=1}^n p*f(p)sum_{d=1}^{lfloorfrac{n}{p} floor}mu(d)sum_{i=1igwedge d|i}^{lfloorfrac{n}{p} floor}sum_{j=1igwedge d|j}^{lfloorfrac{m}{p} floor} ij)

(})

(=sum_{p=1}^n p*f(p)sum_{d=1}^{lfloorfrac{n}{p} floor}mu(d)sum_{i=1}^{lfloorfrac{n}{pd} floor}sum_{j=1}^{lfloorfrac{m}{pd} floor} ijd^2)

(=sum_{p=1}^n p*f(p)sum_{d=1}^{lfloorfrac{n}{p} floor}mu(d)*{d^2}sum({lfloorfrac{n}{pd} floor})sum({lfloorfrac{m}{pd} floor}))

——————————————————————————————————

令(k=pd)，则(p|k)，那么我们可以把sum函数提出并将式子整理得：

(sum_{k=1}^n sum({lfloorfrac{n}{k} floor})sum({lfloorfrac{m}{k} floor})*sum_{p|k} p*f(p)*mu(frac{k}{p})*{({frac{k}{p})}^2})

设(h(k)=sum_{p|k} p*f(p)*mu(frac{k}{p})*{({frac{k}{p})}^2})，则可化简为：

(sum_{k=1}^n sum({lfloorfrac{n}{k} floor})sum({lfloorfrac{m}{k} floor})*h(k))

显然(h(k))之前的那部分可以通过分块来实现(O(sqrt{n}))求解，那么现在关键问题是求(h(k))。

根据各种迷之推导我们得出(h(k))为积性函数，则我们可以用线性筛来求(h(k))。

分情况讨论，当(k)为质数的时候，显然带入解得(h(k)=k(1-k))。

当(k)不为质数的时候，我们可以表示为(k=x*p)（p为质数），如果(x)与(p)互质的话可以直接相乘，否则：

首先对(x)分解质因数为(x={p1^{q1}}{p2^{q2}}{p3^{q3}}...{pk^{qk}})，假设(p=p1)，那么显然：

当(q1>=2)时(f(x*p1)=0)，显然此时(h(x*p1)=0)。

当(q1=1)时(h(x*p1)=h(x/p1)*h(p1^2))(显然(x/p1)和(p1^2)互质)

而(h(p^2)=1*f(1)*mu(p^2)*{p^4}+p*f(p)*mu(1)*{p^2}+{p^2}*f({p^2})*mu(1)*1)

(=0-3p+0=-3p)

( herefore h(x*p1)=-3p1*h(x/p1))

#include<cstdio>
#include<queue>
#include<cctype>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
const int N=4000010;
const int p=(1<<30)-1;
int f[N],su[N],g[N];
bool he[N];
inline int s(int x){
    return x*(x+1)>>1;
}
void Euler(int n){
    int tot=0;
    f[1]=1;    
    for(int i=2;i<=n;i++){    
	if(!he[i]){    
	    su[++tot]=i;
	    f[i]=i*(1-i);
	}
	for(int j=1;j<=tot;j++){
	    int t=i*su[j];
	    if(t>n)break;
	    he[t]=1;
	    if(i%su[j]==0){
		int x=i/su[j];
		if(x%su[j])f[t]=-su[j]*su[j]*su[j]*f[x];
		else f[t]=0;
		break;
	    }
	    else f[t]=f[i]*f[su[j]];
	}
    }
    for(int i=1;i<=n;i++)f[i]+=f[i-1];
    return;
}   
int main(){
    Euler(4000000);
    int t;
    scanf("%d",&t);
    while(t--){
	int n,m,ans=0;
	scanf("%d%d",&n,&m);
	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+=(f[j]-f[i-1])*s(n/i)*s(m/i);
	}
	printf("%d
",ans&p);
    }
    return 0;
}