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}^n\sum_{j=1}^n(n-i+1)(n-j+1)(gcd(i,j)-1)
\]
我们现在只是枚举了斜率为正的情况，而斜率为负数是个镜像问题，方案数与为正相同，所以前面还要乘以2.
这样枚举的复杂度是
\[O(n^2)
\]
的，对于和这个题差不多的一道CQOI的题复杂度是够用的。但是对于本题，
\[n\leqslant 2\cdot10^6
\]
是远远不够的，所以我们必须把算法进一步优化，我这里用的是一种笨方法，强行把这个式子拆成8份。
记：
\[A=(n+1)^2\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)
\]
\[B=\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)
\]
\[C=-(n+1)\sum_{i=1}^n\sum_{j=1}^nigcd(i,j)
\]
\[D=-(n+1)\sum_{i=1}^n\sum_{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}^n\sum_{j=1}^ngcd(i,j)=\sum_{i=1}^n\sum_{j=1}^n\sum_{d|gcd(i,j)}\phi(d)=\sum_{d=1}^n\phi(d)\sum_{i=1}^n[d|i]\sum_{j=1}^n[d|j]=\sum_{d=1}^n\phi(d)\lfloor\frac nd\rfloor\lfloor\frac nd\rfloor
\]
\[\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)=\sum_{d=1}^n\phi(d)\sum_{i=1}^ni[d|i]\sum_{j=1}^nj[d|j]=\sum_{d=1}^nd^2\phi(d)S_1(n/d)S_1(n/d)
\]
\[\sum_{i=1}^n\sum_{j=1}^nigcd(i,j)=\sum_{i=1}^n\sum_{j=1}^njgcd(i,j)=\sum_{d=1}^nd\phi(d)S_1(n/d)
\]
可以直接扫一遍计算，也可以除法分块
\[O(\sqrt n)
\]
计算。
事实上，这题的数据还可以开到更大，例如
\[n\leqslant 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\n", ans);
return 0;
}
