CF#462 div1 D：A Creative Cutout
CF#462 div1 D：A Creative Cutout
题目大意：
原网址戳我！
题目大意：
在网格上任选一个点作为圆中心，然后以其为圆心画\(m\)个圆。
其中第\(k\)个圆的半径为\(\sqrt{k}\),\(m\)个圆的编号依次为\(1\)~\(m\)。
定义一个格点的美妙值\(g(x,y)\)为包含了它的所有圆的编号之和。
定义\(f(n)\)为：当画了\(n\)个圆时，\(f(n) = \sum_{i,j\in R} g(i,j)\)。
现在非常变态的问你一个非常无聊的恶心问题：
给定一个\(m\)(\(m\leq 10^{12}\))，请计算\(ans = \sum_{i = 1}^m f(i)\ mod\ 1000000007\)
样例与样例解释
样例输入输出
\(Input\)：5
\(Output\)：387
样例说明：
样例中，当\(n = 5\)时圆的分布情况如下图所示：
当\(n = 5\)时，如图
有\(5\)个红色的点被圆1、2、3、4、5包含 ， 它们的\(r = \sum_{i = 1}^5 i = 15\)
有\(4\)个橘色的点被圆2、3、4、5包含，它们的\(r = \sum_{i = 2}^5 i = 14\)
有\(4\)个绿色的点被圆4、5包含，它们的\(r = \sum_{i = 4}^5 i = 9\)
有\(8\)个蓝色的点被圆5包含 ， 它们的\(r = 5\)
剩下的灰色的点则不被任意一个圆包含，它们的\(r = 0\)
综上，\(f(5) = 5\times 15 + 4\times 14 + 4\times 9 + 8\times 5 = 207\)
同理可得\(f(1)\)=\(5\)、\(f(2)\)=\(23\)、\(f(3)\)=\(50\)、\(f(4)\)=\(102\)。
所以当\(m = 5\)时 ， \(ans = \sum_{i = 1}^5 f(i) = 5+23+50+102+207 = 387\)
解法与思路
这TM比 \(E\)题 还要难好多好吗？
这里方便叙述我们以圆心为原点建立坐标轴。
对于一个\(f(n)\)，我们可以计算出每一个坐标的贡献：
\(res_n = \sum_{k = x^2 + y^2} ^ n k = \frac{(n+(x^2+y^2))(n-(x^2+y^2)+1)}{2}=\frac{n(n+1) - (x^2+y^2)(x^2+y^2-1)}{2}\)
转换一下得到：\(res = \binom{n+1}{2} - \binom{x^2+y^2}{2}\)
令\(L = x^2 + y^2\)，则每一个\(L\)的答案的贡献为：
\(Res_L = \sum_{k = L}^m (\binom{k+1}{2} - \binom{x^2+y^2}{2}) = \sum_{k = L}^m \binom{k+1}{2} - (m - L +1)\binom{L}{2}\)
然后这里有一个组合公式( 提示 )：
公式：\(\sum_{k = L}^R \binom{k}{w} = \binom{R+1}{w+1} - \binom{L}{w+1}\)
用归纳法证明：
当\(L = R = 1\)时，\(\sum_{k=1}^1\binom{1}{1} = \binom{2}{2} - \binom{1}{2} = 1\)成立。
\(\sum_{k = L}^{R+1}\binom{k}{w} = \binom{R+1}{w}+\sum_{k=L}^R\binom{k}{w} = [\binom{R+1}{w+1}+ \binom{R+1}{w}] - \binom{L}{w+1}\)
由组合数的计算公式可得\(\binom{R+1}{w+1}+ \binom{R+1}{w} = \binom{R+2}{w+1}\)
所以\(\sum_{k = L}^{R+1} \binom{k}{w} = \binom{R+2}{w+1} - \binom{L}{w+1}\)依旧成立。
同理可以证明\(L\)变化到\(L+1\)此结论任符合。
综上：\(\sum_{k = L}^R \binom{k}{w} = \binom{R+1}{w+1} - \binom{L}{w+1}\)
我们套用这个公式：
\(Res_L = \binom{m+2}{3} - \binom{L+1}{3} - (m-L+1)\binom{L}{2}\)
\(Res_L = \frac{1}{6}(\frac{(m+2)!}{(m-1)!} - \frac{(L+1)!}{(L-2)!} - (m-L+1)\frac{3L!}{(L-2)!})\)
暴力化简一顿之后得到：
嗯，一个关于\(L\)的多项式。
考虑一下枚举\(x^2\) 与 \(y^2\)的话：
令\(L = x^2 + y^2\)的话，带入原式中可以得到：
\(Res_L = \frac{1}{6}(\)
\(\ \ \ \ \ 2*y^6+\)
\(\ \ \ \ (6x^2-(3m+6))*y^4+\)
\(\ \ \ \ (6x^4-2(3m+6)x^2+(3m+4))*y^2+\)
\(\ \ \ \ (2x^6-(3m+6)x^4+(3m+4)x^2+(m)(m+1)(m+2))*y^0\)
\()\)
等于说，我们如果枚举了\(x\)，那么上面除\(y\)以外其它的都是系数(常数)。
我们把它命名为一个关于\(y\)的函数\(S_x(y)\)，那么答案可以写为：
所以说，我们只需要枚举\(x\)，然后计算\(\sum_{y\in Z}^{y^2 \leq m - x^2} S_x(y)\)即可。
而\(S_x(y)\)中，我们只需要计算\(y^2\)、\(y^4\)、\(y^6\)即可。
然后对于上面这三项，我们单独算答案，可以直接套数学公式：
\(\sum_{i=1}^{N}i^{2}=\frac{N(N+1)(2N+1)}{6}\)
\(\sum_{i=1}^{N}i^{4}=\frac{N(N+1)(2N+1)(3N^2+3N-1)}{30}\)
\(\sum_{i=1}^{N}i^{6}=\frac{N(N+1)(2N+1)(3N^4+6N^3-3N+1)}{42}\)
求\(\sum_{i=1}^N i^{2t}\)的\(O(1)\)公式，用\((r+1)^{2t+1} - r^{2t+1}\)变形即可得到。
枚举\(x\)只需要枚举\(\sqrt{m}\)次，所以总的复杂度为\(O(\sqrt{m})\)。
实现代码：
注意不要暴\(int64\)了!!
枚举\(x\)直接从\(-\sqrt{m}\)枚举到\(\sqrt{m}\)即可，省去讨论的麻烦。
#include<bits/stdc++.h>
#define ll long long
#define MOD 1000000007
using namespace std;
ll mod(ll e){ e %= MOD; if(e < 0)e += MOD; return e; }
ll inv6 , inv30 , inv42;
ll f[5] , n , m , x , y , x2 , x4 , x6 , ans;
ll Pow(ll T , ll js , ll S){
while(js){ if(js&1)S = mod(S * T);
T = mod(T * T); js >>= 1; } return S;
}
ll sum(ll i , ll k){
if(k == 0)return mod( 2*i + 1 );
ll s0 = mod( 2 * mod( mod( i * (i+1) ) * mod(2*i + 1) ) );
if(k == 2)return mod(inv6 * s0);
ll t0 , t1 , t2 , t3 , t4;
t0 = 1; t1 = i; t2 = mod(i * i); t3 = mod(t1 * t2); t4 = mod(t2 * t2);
if(k == 4)return mod(inv30 * mod(s0 * mod( mod(3 * t2) + mod(3 * t1) - t0 )));
if(k == 6)return mod(inv42 * mod(s0 * mod(mod(3*t4) + mod(6*t3) - mod(3*t1) + t0) ) );
}
int main(){
ans = 0; cin >> m; n = sqrt(m);
f[0] = mod( mod(mod(m) * mod(m+1)) * mod(m+2) ); //注意这里别爆int64了！
f[1] = mod( 3*m + 4 );
f[2] = mod( -1 * (3*m + 6) );
f[3] = 2;
inv6 = Pow(6 , MOD - 2 , 1);
inv30 = Pow(30 , MOD - 2 ,1);
inv42 = Pow(42 , MOD - 2 , 1);
for(ll x = -n; x <= n; x ++){
y = sqrt(m - x * x);
x2 = mod(x * x); x4 = mod(x2 * x2); x6 = mod(x2 * x4);
ans = mod( ans + mod(f[3] * sum(y , 6)) );
ans = mod( ans + mod( mod( mod(6 * x2) + f[2] ) * sum(y , 4) ) );
ans = mod( ans + mod( mod(mod(6 * x4) + mod(mod(2 * f[2])*x2) + f[1]) * sum(y , 2)) );
ans = mod( ans + mod(mod(f[3]*x6) + mod(f[2]*x4) + mod(f[1]*x2) + f[0]) * sum(y , 0) );
ans = mod( ans );
}
ans = mod( ans * inv6 ); cout << ans; return 0;
}
