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$個紅色的點被圓一、二、三、四、5包含 ， 它們的$r = \sum_{i = 1}^5 i = 15$ 有$4$個橘色的點被圓二、三、四、5包含，它們的$r = \sum_{i = 2}^5 i = 14$ 有$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)!})$ 暴力化簡一頓以後獲得： $$Res_L = \frac{1}{6}[\ 2L^3-3(m+2)L^2+(3m+4)L+m(m+1)(m+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)$，那麼答案能夠寫爲： $$Ans = \frac{1}{6}\sum_{x \in Z }^{x^2 \leq m}\ \sum_{y\in Z}^{y^2 \leq m - x^2} 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;
}