P3166 [CQOI2014]数三角形

题目

P3166 [CQOI2014]数三角形

分析

数三点不共线显然不如数三点共线（其实是因为在这个之前做了一道不共线的 BZOJ3518点组计数。）

首先除去边角的三点共线，因为很好算，我们就可以只考虑斜着的三点共线。

直接枚举两个点代价太大，而枚举一个点信息又太少，于是考虑枚举第一个点和第三个点的横纵坐标之差，考虑有这样的差的对数有多少个呢？ \((n-i)(m-j)\) 个 。

这样第二个点的个数就可以确定了是 \(\gcd(i,j)-1\) ，直接乘起来然后求和即可。

现在问题转化成求 \large \sum\limits_{i=1}^{n\sum\limits_{j=1}}m(gcd(i,j)-1)(n-i)(m-j) 的值。

\large \sum\limits_{i=1}^{n\sum\limits_{j=1}}m(gcd(i,j)-1)(n-i)(m-j)

\large =\sum\limits_{i=1}^{n\sum\limits_{j=1}}m(gcd(i,j))(n-i)(m-j)-\sum\limits_{i=1}^{n{\sum\limits_{j=1}}m{ij}}

\large =\sum\limits_{d=1}^{{min(n,m)}{\varphi{(d)}\sum\limits_{i=1}}\rfloor}{(n-di)\sum\limits_{j=1}^{\lfloor \dfrac{m}{d}\rfloor}(m-dj)}}-\sum\limits_{i=1}^{{n-1}{i\sum\limits_{j=1}}{j}}

因为 n,m 很小，所以可以直接暴力计算。（注意中间两个求和柿子就是等差数列，我当时竟然没看出来。。）

最后使用减法原理即可。

代码

#include<bits/stdc++.h>
using namespace std;
template <typename T>
inline void read(T &x){
    x=0;char ch=getchar();bool f=false;
    while(!isdigit(ch)){if(ch=='-'){f=true;}ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return ;
}
template <typename T>
inline void write(T x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
    return ;
}
#define ll long long
#define ull unsigned long long
const int N=1e6+5,M=1e6+5,MOD=1e9+7;
ll n,m;
ll cnt,prime[N],f[N],low[N];
bool vis[N];
inline void GetPrimes(int d){
	vis[1]=true;f[1]=low[1]=1;//对1进行定义 
	for(ll i=2;i<=d;i++){
		if(!vis[i]) prime[++cnt]=low[i]=i,f[i]=i-1;//对质数进行定义 
		for(ll j=1;j<=cnt&&i*prime[j]<=d;j++){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
                low[i*prime[j]]=low[i]*prime[j];
                if (low[i]==i) f[i*prime[j]]=f[i]*prime[j];//对质数的若干次幂进行定义(一般由f[i]递推);
                else f[i*prime[j]]=f[i/low[i]]*f[low[i]*prime[j]];
                break;
            }
            low[i*prime[j]]=prime[j];
			f[i*prime[j]]=f[i]*f[prime[j]];
		}
	}
	return ;
}
int main(){
	GetPrimes(1001);
	read(n),read(m);n++,m++;
	ll d=min(n,m);ll Ans=0;
	for(ll i=1;i<=d;i++) Ans=(Ans+f[i]*((n-i+n%i)*((int)floor(n/i))/2)*((m-i+m%i)*((int)floor(m/i))/2));
	Ans-=n*(n-1)/2*m*(m-1)/2;
	Ans=(Ans+Ans);
	Ans=(Ans+n*(m*(m-1)*(m-2)/6)+m*(n*(n-1)*(n-2)/6));
	Ans=n*m*(n*m-1)*(n*m-2)/6-Ans;
	write(Ans);
	return 0;
}