校内互测 光 （数论+容斥原理）

光

【问题描述】
天猫有一个长方形盒子，长宽分别为A，B。这个长方形盒子的内壁全部是镜面。天猫在这个盒子的左下方放了一个激光灯。这个灯可以照向盒子内的任意角度。现在天猫想要打开这个激光灯，但是他想让光线按照如下规则照射：
1.这束光必须恰好打到盒子边缘反射D次，并且不能碰到任意一个角落（除了出发点以及结束点）。
2.这束光必须到达盒子右上角，并且结束反射。天猫想要知道，所有合法的光线路线的长度平方和是多少。作为一个
资深OIer，你应该知道输出要对10^9+7取模。
【输入文件】
一行三个数，表示A，B，D。
【输出文件】
一行表示路线长度的平方和。
【样例输入】
[sample 1]
3 4 0
[sample 2]
3 3 2
【样例输出】
[sample 1]
25
[sample 2]
180
【数据范围】
       对于20%的数据，D<=2
 对于40%的数据，D<=10
 对于70%的数据，D<=10^6
 对于100%的数据，D<=10^9, A,B<=10^6

题解：数论+容斥原理

注意这个题的意思是在平面内不断的发生镜面反射，让你寻找从左下角到右上角的路径。我们可以倒着每次找出对应出入射点的虚像，然后连线，不断把得到的直线关于镜面的边缘对称过去，最终一定会得到一条过左下角的直线。

那么问题其实可以转换成有一个d*d的网格，我们要在网格中选择一个点(x,y),要满足x+y-2=d，其实(0,0)到(x,y)的连线经过的格子的边界就是每次入射的位置，那么我们要求在路径中不能经过任何角落，所以我们需要保证这条直线不经过任何的整点，也就是gcd(x,y)=1.

那么答案就是下面式子的二分之一，因为每个点都计算了两次，因为(x^2+y^2)*(a^2+b^2)得到的实际上是个四项式，也就是(x,y)与(y,x)的答案。

那么这个问题解决的关键就是如何求1到n内与n互质的数的平方和（发现上面的式子是之所以算了两次，是因为每个合法的i^2都算了两次，其实我们可以直接计算1到n内与n互质的数的平方和，这个一定是两两对应形成点对）。

我们设n=d+2，那么其实可以用容斥原理来解决这个问题，我们要求的数是满足gcd(x,n)=1的数，也就是与n相同的质因子的个数是0的时候。ans=至少有0个与n相同的质因子-至少有1个与n相同的质因子+至少有2个与n相同的质因子....

至少有0个的其实就是1到n的平方和，有公式n(n+1)(2*n+1)/6

至少有1个，那么x必然是p(质因子)的倍数，也就是sigma(i=1..(n/p)) i*p ,这些数的平方和其实就是p^2*sigma(i=1..n/p)i^2 

然后写dfs求解就可以了。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<algorithm>
#define N 100000
#define LL long long 
#define p 1000000007
using namespace std;
LL prime[N+10],pd[N+10],num[N],cnt;
LL a,b,d,ans,n;
LL pow(LL x)
{
	return x*x%p;
}
LL mul(LL num,LL x)
{
	LL ans=1; LL base=num%p;
	while (x){
	   if (x&1) ans=ans*base%p;
	   x>>=1;
	   base=base*base%p;
	}
	return ans;
}
LL sqrsum(LL x)
{
	LL t=x*(x+1)%p*(2*x+1)%p*mul(6,p-2)%p;
	//cout<<x<<" "<<t<<endl;
	return t%p;
}
void calc()
{
	for (int i=2;i<=N;i++) {
		if (!pd[i]) {
		   prime[++prime[0]]=i;
		}
		for (int j=1;j<=prime[0];j++){
			if (prime[j]*i>N) break;
			pd[prime[j]*i]=1;
			if (i%prime[j]==0) break;
		}
	}
}
void solve(LL d)
{
	for (int i=1;prime[i]*prime[i]<=d;i++)
	if (d%prime[i]==0){
		num[++cnt]=prime[i];
		while (d%prime[i]==0) d/=prime[i];
	}
	if (d>1) num[++cnt]=d;
}
LL getf(LL i)
{
	LL size=n/i;
	LL ans=sqrsum(size)*pow(i); 
	return ans%p;
}
void dfs(LL tot,int x,int now)
{
	if (now>cnt){
		//if (x==0) return;
		if (x%2)  ans-=getf(tot); 
		else ans+=getf(tot);
		ans=(ans%p+p)%p;
		//cout<<ans<<endl;
		return;
    }
	dfs(tot*num[now]%p,x+1,now+1);
	dfs(tot,x,now+1);
}
int main()
{
	freopen("light.in","r",stdin);
	freopen("light.out","w",stdout);
	calc();
	scanf("%I64d%I64d%I64d",&a,&b,&d);
	if (d&1) {
		printf("0\n");
		return 0;
	}
	LL dis=(pow(a)+pow(b))%p; 
	n=d+2; 
	solve(n); 
	//for (int i=1;i<=cnt;i++) cout<<num[i]<<" ";
	//cout<<endl;
	dfs(1,0,1); //cout<<ans<<endl;
	printf("%I64d\n",dis*ans%p);
}