EOJ 数三角形
3618. 數三角形
題面統計數據討論
單點時限: 2.0 sec
內存限制: 256 MB
求在 n×n 矩形內(包括邊界)整點三角形有多少個。由于答案可能很大,對于 109+7 取模。
輸入格式
第一行一個整數 n (1≤n≤2?106) 表示矩形的寬與高。
輸出格式
一個整數 x。
樣例
input
2
output
76
input
5
output
6768
提示
做不出來可以試試打表找規律?
題解:這個題貌似只有4個人在eoj做出來。。。看著也挺難的,光氣場就難住了。。。
但其實,這道題并不難,就是有一些繁瑣罷了(甚至可以說連繁瑣都算不上)。
分析一下,首先這是一道推數學公式的題,數格點三角形。我們很容易想到在格點上能構成三角形的條件是三個頂點不共線,其余的情況都可以構成三角形,原因很顯然,在平面上畫三個點能構成三角形的條件也是如此。我們發現正著做不好做,于是我們就反著做,考慮三點共線的情況:只有三種情況,三點豎著,三點橫著,三點斜著,橫豎很好用組合數解決,于是現在就只有斜著的了。
考慮斜著的兩個點,中間不會夾著其他點的充要條件是
[gcd(x,y)=1
]
其中,x,y分別為兩點間橫縱坐標差。
當兩個點夾著中間點的情況,點數恰好為
[gcd(x,y)-1
]
所以這樣我們初步可以用四個計數變量去算這個式子(因為一個點就有兩個維度),但是顯然不行。我們發現有很多情況是相似的,點的位置只是平移了一段,所以我們不妨直接將靠下的點設為(0,0),然后最上的點可以用兩個計數變量枚舉,可以設為(i,j),然后兩點形成的矩形可以嵌入到正方形中,共有
[(n-i+1)(n-j+1)
]
種方案。
嵌入的方案中再去選中間點,一共有
[(n-i+1)(n-j+1)(gcd(i,j)-1)
]
種方案。
合起來就是
[sum_{i=1}^nsum_{j=1}^n(n-i+1)(n-j+1)(gcd(i,j)-1)
]
我們現在只是枚舉了斜率為正的情況,而斜率為負數是個鏡像問題,方案數與為正相同,所以前面還要乘以2.
這樣枚舉的復雜度是
[O(n^2)
]
的,對于和這個題差不多的一道CQOI的題復雜度是夠用的。但是對于本題,
[nleqslant 2cdot10^6
]
是遠遠不夠的,所以我們必須把算法進一步優化,我這里用的是一種笨方法,強行把這個式子拆成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都是常數,可以直接算,對于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)
]
可以直接掃一遍計算,也可以除法分塊
[O(sqrt n)
]
計算。
事實上,這題的數據還可以開到更大,例如
[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;
}
總結
- 上一篇: 颧骨内推手术危险吗
- 下一篇: 中国实木家具十大品牌