3618. 數(shù)三角形

題面統(tǒng)計(jì)數(shù)據(jù)討論

單點(diǎn)時(shí)限: 2.0 sec

內(nèi)存限制: 256 MB

求在 n×n 矩形內(nèi)（包括邊界）整點(diǎn)三角形有多少個(gè)。由于答案可能很大，對(duì)于 109+7 取模。

輸入格式

第一行一個(gè)整數(shù) n (1≤n≤2?106) 表示矩形的寬與高。

輸出格式

一個(gè)整數(shù) x。

樣例

input

2

output

76

input

5

output

6768

提示

做不出來(lái)可以試試打表找規(guī)律？

題解：這個(gè)題貌似只有4個(gè)人在eoj做出來(lái)。。?？粗餐﹄y的，光氣場(chǎng)就難住了。。。

但其實(shí)，這道題并不難，就是有一些繁瑣罷了（甚至可以說(shuō)連繁瑣都算不上）。

分析一下，首先這是一道推數(shù)學(xué)公式的題，數(shù)格點(diǎn)三角形。我們很容易想到在格點(diǎn)上能構(gòu)成三角形的條件是三個(gè)頂點(diǎn)不共線，其余的情況都可以構(gòu)成三角形，原因很顯然，在平面上畫(huà)三個(gè)點(diǎn)能構(gòu)成三角形的條件也是如此。我們發(fā)現(xiàn)正著做不好做，于是我們就反著做，考慮三點(diǎn)共線的情況：只有三種情況，三點(diǎn)豎著，三點(diǎn)橫著，三點(diǎn)斜著，橫豎很好用組合數(shù)解決，于是現(xiàn)在就只有斜著的了。

考慮斜著的兩個(gè)點(diǎn)，中間不會(huì)夾著其他點(diǎn)的充要條件是

[gcd(x,y)=1
]

其中，x,y分別為兩點(diǎn)間橫縱坐標(biāo)差。

當(dāng)兩個(gè)點(diǎn)夾著中間點(diǎn)的情況，點(diǎn)數(shù)恰好為

[gcd(x,y)-1
]

所以這樣我們初步可以用四個(gè)計(jì)數(shù)變量去算這個(gè)式子（因?yàn)橐粋€(gè)點(diǎn)就有兩個(gè)維度），但是顯然不行。我們發(fā)現(xiàn)有很多情況是相似的，點(diǎn)的位置只是平移了一段，所以我們不妨直接將靠下的點(diǎn)設(shè)為(0,0)，然后最上的點(diǎn)可以用兩個(gè)計(jì)數(shù)變量枚舉，可以設(shè)為(i,j)，然后兩點(diǎn)形成的矩形可以嵌入到正方形中，共有

[(n-i+1)(n-j+1)
]

種方案。

嵌入的方案中再去選中間點(diǎn)，一共有

[(n-i+1)(n-j+1)(gcd(i,j)-1)
]

種方案。

合起來(lái)就是

[sum_{i=1}^nsum_{j=1}^n(n-i+1)(n-j+1)(gcd(i,j)-1)
]

我們現(xiàn)在只是枚舉了斜率為正的情況，而斜率為負(fù)數(shù)是個(gè)鏡像問(wèn)題，方案數(shù)與為正相同，所以前面還要乘以2.

這樣枚舉的復(fù)雜度是

[O(n^2)
]

的，對(duì)于和這個(gè)題差不多的一道CQOI的題復(fù)雜度是夠用的。但是對(duì)于本題，

[nleqslant 2cdot10^6
]

是遠(yuǎn)遠(yuǎn)不夠的，所以我們必須把算法進(jìn)一步優(yōu)化，我這里用的是一種笨方法，強(qiáng)行把這個(gè)式子拆成8份。

記：

[A=(n+1)^2sum_{i=1}^nsum_{j=1}^ngcd(i,j)
]

[B=sum_{i=1}^nsum_{j=1}^nijgcd(i,j)
]

[C=-(n+1)sum_{i=1}^nsum_{j=1}^nigcd(i,j)
]

[D=-(n+1)sum_{i=1}^nsum_{j=1}^njgcd(i,j)
]

[E=n^2(n+1)^2
]

[F=S_3(n),S_3(n)=sum_{k=1}^nk^3
]

[G=-n(n+1)S_1(n),S_1(n)=sum_{k=1}^nk
]

[H=-n(n+1)S_1(n),S_1(n)=sum_{k=1}^nk
]

答案就是

[A+B+C+D-E-F-G-H
]

其中E,F,G,H都是常數(shù)，可以直接算，對(duì)于A,B,C,D

[sum_{i=1}^nsum_{j=1}^ngcd(i,j)=sum_{i=1}^nsum_{j=1}^nsum_{d|gcd(i,j)}phi(d)=sum_{d=1}^nphi(d)sum_{i=1}^n[d|i]sum_{j=1}^n[d|j]=sum_{d=1}^nphi(d)lfloorfrac ndfloorlfloorfrac ndfloor
]

[sum_{i=1}^nsum_{j=1}^nijgcd(i,j)=sum_{d=1}^nphi(d)sum_{i=1}^ni[d|i]sum_{j=1}^nj[d|j]=sum_{d=1}^nd^2phi(d)S_1(n/d)S_1(n/d)
]

[sum_{i=1}^nsum_{j=1}^nigcd(i,j)=sum_{i=1}^nsum_{j=1}^njgcd(i,j)=sum_{d=1}^ndphi(d)S_1(n/d)
]

可以直接掃一遍計(jì)算，也可以除法分塊

[O(sqrt n)
]

計(jì)算。

事實(shí)上，這題的數(shù)據(jù)還可以開(kāi)到更大，例如

[nleqslant 10^9
]

然后用杜教篩求解。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;//simplify long long
typedef unsigned long long ull;
#define inf 2147483647
#define pi 3.14159265358979
#define rep(i, l, r) for(int i = l; i <= r; i ++)
#define lop(i, r, l) for(int i = r; i >= l; i --)
#define step(i, l, r, __step) for(int i = l; i <= r; i += __step)
#define revp(i, r, l, __step) for(int i = r; i >= l; i -= __step)
#define regsiter reg
#define regsiter int RI
#define regsiter long long RL
inline int read()//fast read
{
	int ret = 0, sgn = 1;
	char chr = getchar();
	while(chr < '0' || chr > '9')
	{if(chr == '-') sgn = -1; chr = getchar();}
	while(chr >= '0' && chr <= '9')
	{ret = ret * 10 + chr - '0'; chr = getchar();}
	return ret * sgn;
}
const int N = 2e6 + 50;
const ll mod = 1e9 + 7;
ll inv6;
int prime[N], cnt, phi[N];
ll sumphi[N], sum2[N], sum1[N];
bool vis[N];
void shai(int x)
{
	phi[1] = 1;
	rep(i, 2, x)
	{
		if(!vis[i])
		{
			prime[++ cnt] = i;
			phi[i] = i - 1;
		}
		rep(j, 1, cnt)
		{
			if(i * prime[j] > x) break;
			vis[i * prime[j]] = 1;
			if(!(i % prime[j]))
			{
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
			phi[i * prime[j]] = phi[i] * (prime[j] - 1);
		}
	}
	rep(i, 1, x)
	{
		sumphi[i] = sumphi[i - 1] + phi[i], sumphi[i] %= mod;
		sum2[i] = sum2[i - 1] + 1ll * i * i % mod * phi[i] % mod, sum2[i] %= mod;
		sum1[i] = sum1[i - 1] + 1ll * i * phi[i] % mod, sum1[i] %= mod;
	}
}
ll ksc(ll x, ll y)
{
	ll ret = 0;
	while(y)
	{
		if(y & 1) ret += x, ret %= mod;
		x += x, x %= mod;
		y >>= 1;
	}
	return ret;
}
ll ksm(ll x, ll y)
{
	ll ret = 1;
	while(y)
	{
		if(y & 1) ret *= x, ret %= mod;
		x *= x, x %= mod;
		y >>= 1;
	}
	return ret;
}
ll s1(ll x)
{
	return (x * (x + 1) / 2) % mod;
}
ll s2(ll x)
{
	return x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
}
ll s3(ll x)
{
	return (x * (x + 1) / 2) % mod * ((x * (x + 1) / 2) % mod) % mod;
}
ll c3(ll x)
{
	if(x < 3) return 0;
//	x %= mod;
	return x % mod * (x - 1) % mod * (x - 2) % mod * inv6 % mod;
}
int main()
{
	ll n;
	inv6 = ksm(6, mod - 2);
	shai(N - 10);
	scanf("%lld", &n);
	ll sumA = 0, sumB = 0, sumC = 0;
	for(int l = 1, r; l <= n; l = r + 1)
	{
		r = min(n / (n / l), n);
		sumA += (sumphi[r] - sumphi[l - 1] + mod) % mod * 1ll * (n / l) % mod * 
		(n / l) % mod;
		sumA %= mod;
	}
	sumA *= 1ll * (n + 1) % mod * (n + 1) % mod, sumA %= mod;
	for(int l = 1, r; l <= n; l = r + 1)
	{
		r = min(n / (n / l), n);
		sumB += (sum2[r] - sum2[l - 1] + mod) % mod * s3(n / l) % mod;
		sumB %= mod;
	}
//	cout << sumB << endl;
	for(int l = 1, r; l <= n; l = r + 1)
	{
		r = min(n / (n / l), n);
		sumC += (sum1[r] - sum1[l - 1] + mod) % mod * s1(n / l) % mod * (n / l) % mod;
		sumC %= mod;
	}
	sumC = ksc(sumC, n + 1);
	sumC *= -2;
//	sumC *= -2ll * (n + 1);
	sumC = (sumC % mod + mod) % mod;
//	cout << sumC << endl;
	ll sumE = 1ll * n * n % mod * (n + 1) % mod * (n + 1) % mod;
	sumE %= mod;
	ll sumF = s3(n);
	ll sumG = -sumE;
//	ll sumG = -2ll * n * (n + 1) % mod * s1(n);
	sumG = (sumG % mod + mod) % mod;
	ll SS = (sumA + sumB) % mod;
	SS += sumC, SS %= mod;
	ll TT = (sumE + sumF) % mod;
	TT += sumG, TT %= mod;
//	SS -= sumE, (SS += mod) %= mod;
//	SS -= sumF, (SS += mod) %= mod;
//	SS -= sumG, (SS += mod) %= mod;
	SS -= TT, SS = (SS % mod + mod) % mod;
//	cout << SS << endl;
	SS *= 2, SS %= mod;
	ll m = n + 1;
	SS += 2ll * ksc(m, c3(m)) % mod, SS %= mod;
	SS = (SS % mod + mod) % mod;
	ll P = m * m % mod;
	ll Q = -SS;
	Q = (Q % mod + mod) % mod;
	ll ans = c3(P) - SS;
//	ll ans = ksc(-SS, c3(P));
	ans = (ans % mod + mod) % mod;
	printf("%lld
", ans);
	return 0;
}

總結(jié)

以上是生活随笔為你收集整理的EOJ 数三角形的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。