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\)

解法與思路

這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}\)

公式：\(\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;
}

CF#462 div1 D：A Creative Cutout