【也许是真的最后一篇】炫酷式子题

Posted on:2023年11月9日 at 08:21

本文是一篇连载文章，系列“题解”，同样系列的文章还有这些

这个东西放在做题笔记里完全不能彰显它的毒瘤，我就单独拿出来了。

题意

有 n(n≤1014)n(n \le 10^{14}) 个圆心为原点，半径为 1,2,3,⋯ ,n1, \sqrt{2}, \sqrt{3}, \cdots, \sqrt{n} 的圆，定义 f(i)f(i) 为求出每个点被 ii 个圆覆盖的编号和，你需要求出 ∑i∈[1,n]f(i)\sum_{i \in [1, n]}f(i)（答案对 109+710^9+7 取模）。

样例：5 | 387, 88 | 16605192, 888888 | 202495041, 88888888888 | 143398458

解法

n≤1014n \le 10^{14} 给了我们提示，发现最远被覆盖的点的坐标绝对值不会超过 10710^7，所以我们可以考虑枚举一维，然后 O(1)O(1) 地算出这条线上的答案。

写之前先作一下定义，gk(i)=∑j∈[0,i]jkg_{k}(i) = \sum_{j \in [0, i]}j^k，也就是 kk 次方前缀和。

首先有一个显然的对称性，也就是一个圆可以分成四个象限来算，最后零点的单独算。

（下文仅考虑一象限和 x 轴的正半轴上）

在第一象限和 x 坐标轴的正半轴上，我们枚举 xx 坐标，然后可以得出在这个点上最高的坐标能取到多少。

假定 (x,y)(x, y) 可以被取到，总共有 nn 个圆，则它到原点的距离为 x2+y2\sqrt{x^2+y^2}，会被 x2+y2,x2+y2+1,x2+y2+2,⋯ ,nx^2+y^2, x^2+y^2+1, x^2+y^2+2, \cdots, n 号圆覆盖，然后因为我们要求 ff 的前缀和，发现若被 11 号点覆盖，会重复计算 nn 次，若被 22 号点覆盖，会重复计算 n−1n-1 次，依此类推。

所以我们要求的是 ∑k∈[x2+y2,n]k(n+1−k)\sum_{k \in [x^2+y^2, n]}k(n+1-k)。

在零点上，会被所有圆覆盖，所以答案应额外加上 ∑k∈[1,n]k(n+1−k)=n(n+1)22−g2(n)\sum_{k \in [1, n]}k(n+1-k) = \frac{n(n+1)^2}{2} -g_2(n)。

然后就开始推式子了。

如果公式有误的话，下面有每一行公式的代码实现，可以先看看代码（都是过了样例的，只是复杂度问题），然后联系我纠正。

可能有点啰嗦，但希望大家看懂，还有就是中间的代码可能得多取模几下，后面调了些东西。

这道题搞得我再也不想推式子了，不过马上就退役了。

n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]∑k∈[x2+y2,n]k(n+1−k)=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]∑k∈[x2+y2,n](k(n+1)−k2)=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]((n+1)∑k∈[x2+y2,n]k−∑k∈[x2+y2,n]k2)=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]((n+1)(x2+y2+n)(n−x2−y2+1)2−∑k∈[x2+y2,n]k2)=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]((n+1)(x2+y2+n)(n−x2−y2+1)2+∑k∈[1,x2+y2)k2−∑k∈[1,n]k2)=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]((n+1)(x2+y2+n)(n−x2−y2+1)2+g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]((n+1)(x2+y2−x4−y4−2x2y2+n(n+1))2)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋n(n+1)+∑y∈[0,⌊n−x2⌋](x2+y2−x4−y4−2x2y2))2+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+∑y∈[0,⌊n−x2⌋](y2−y4)−2x2∑y∈[0,⌊n−x2⌋]y2)2+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−∑y∈[0,⌊n−x2⌋]y4−2x2g2(⌊n−x2⌋))2+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](g2(x2+y2−1)−g2(n))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2+4∑x∈[1,⌊n⌋](−⌊n−x2+1⌋g2(n)∑y∈[0,⌊n−x2⌋]g2(x2+y2−1))=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2−4∑x∈[1,⌊n⌋]⌊n−x2+1⌋g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](x2+y2−1)(x2+y2)(2x2+2y2−1)6=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2−4∑x∈[1,⌊n⌋]⌊n−x2+1⌋g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋](x4+y4+2x2y2)(2x2+2y2−1)6=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2−4∑x∈[1,⌊n⌋]⌊n−x2+1⌋g2(n)+4∑x∈[1,⌊n⌋]∑y∈[0,⌊n−x2⌋]2x6+2y6+6x2y4+6x4y2+x2+y2−6x2y2−3x4−3y46=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2−4∑x∈[1,⌊n⌋]⌊n−x2+1⌋g2(n)+4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(2x6+x2)+∑y∈[0,⌊n−x2⌋]2y6+6x2y4+6x4y2+y2−6x2y2−3x4−3y46)=n(n+1)22−g2(n)+(n+1)4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(x2−x4+n(n+1))+g2(⌊n−x2⌋)−g4(⌊n−x2⌋)−2x2g2(⌊n−x2⌋))2−4∑x∈[1,⌊n⌋]⌊n−x2+1⌋g2(n)+4∑x∈[1,⌊n⌋](⌊n−x2+1⌋(2x6+x2−3x4)+−6x2g2(⌊n−x2⌋)+6x2g4(⌊n−x2⌋)+2g6(⌊n−x2⌋)+6x4g2(⌊n−x2⌋)+g2(⌊n−x2⌋)−3g4(⌊n−x2⌋)6)\begin{aligned} % Iteration 1 &\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \sum_{k \in [x^2+y^2, n]} k \left( n+1-k \right) \\ % Iteration 2 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \sum_{k \in [x^2+y^2, n]} \left( k \left( n+1 \right) - k^2 \right) \\ % Iteration 3 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( \left( n+1 \right) \sum_{k \in [x^2+y^2, n]} k - \sum_{k \in [x^2+y^2, n]} k^2 \right) \\ % Iteration 4 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( \frac{\left( n+1 \right) \left( x^2+y^2+n \right) \left( n-x^2-y^2+1 \right)}{2} - \sum_{k \in [x^2+y^2, n]} k^2 \right) \\ % Iteration 5 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( \frac{\left( n+1 \right) \left( x^2+y^2+n \right) \left( n-x^2-y^2+1 \right)}{2} + \sum_{k \in [1, x^2+y^2)} k^2 - \sum_{k \in [1, n]} k^2 \right) \\ % Iteration 6 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( \frac{\left( n+1 \right) \left( x^2+y^2+n \right) \left( n-x^2-y^2+1 \right)}{2} + g_2(x^2+y^2-1) - g_2(n) \right) \\ % Iteration 7 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( \frac{\left( n+1 \right) \left( x^2 + y^2 - x^4 - y^4 - 2x^2y^2 + n(n+1) \right)}{2} \right) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( g_2(x^2+y^2-1)- g_2(n)\right) \\ % Iteration 8 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor n(n+1) + \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( x^2 + y^2 - x^4 - y^4 - 2x^2y^2 \right) \right)}{2} + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( g_2(x^2+y^2-1)- g_2(n)\right) \\ % Iteration 9 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( y^2 - y^4 \right) -2x^2 \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} y^2 \right)}{2} + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( g_2(x^2+y^2-1)- g_2(n)\right) \\ % Iteration 10 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} y^4 -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( g_2(x^2+y^2-1)- g_2(n)\right) \\ % Iteration 11 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]}\sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \left( g_2(x^2+y^2-1)- g_2(n)\right) \\ % Iteration 12 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \left( -\lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} g_2(x^2+y^2-1) \right) \\ % Iteration 13 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} -4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \frac{\left(x^2+y^2-1\right) \left(x^2+y^2\right) \left(2x^2+2y^2-1\right)}{6} \\ % Iteration 14 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} -4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \frac{\left( x^4+y^4+2x^2y^2 \right) \left( 2x^2+2y^2-1 \right)}{6} \\ % Iteration 15 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} -4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \frac{2x^6 + 2y^6 + 6x^2y^4 + 6x^4y^2 + x^2 + y^2 -6x^2y^2 - 3x^4 - 3y^4}{6} \\ % Iteration 16 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} -4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \left( \lfloor \sqrt{n-x^2} + 1 \rfloor (2x^6+x^2) + \sum_{y\in[0, \lfloor\sqrt{n-x^2}\rfloor]} \frac{2y^6 + 6x^2y^4 + 6x^4y^2 + y^2 -6x^2y^2 - 3x^4 - 3y^4}{6} \right) \\ % Iteration 17 =&\frac{n \left( n+1 \right) ^2}{2} - g_2 \left( n \right) + \left( n+1 \right) 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \frac{\left( \lfloor \sqrt{n-x^2}+1 \rfloor \left( x^2 - x^4 + n(n+1) \right) + g_2\left(\lfloor \sqrt{n - x^2} \rfloor \right) - g_4 \left( \lfloor\sqrt{n-x^2}\rfloor \right) -2x^2 g_2\left(\lfloor\sqrt{n-x^2}\rfloor \right) \right)}{2} -4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \lfloor \sqrt{n-x^2}+1 \rfloor g_2(n) + 4\sum_{x \in [1,\lfloor\sqrt{n}\rfloor]} \left( \lfloor \sqrt{n-x^2} + 1 \rfloor (2x^6+x^2-3x^4) + \frac{-6x^2g_2 \left(\lfloor \sqrt{n-x^2} \rfloor \right) + 6x^2g_4 \left(\lfloor \sqrt{n-x^2} \rfloor \right) + 2g_6 \left(\lfloor \sqrt{n-x^2} \rfloor \right) + 6x^4g_2 \left(\lfloor \sqrt{n-x^2} \rfloor \right) + g_2\left(\lfloor \sqrt{n-x^2} \rfloor \right) - 3g_4\left(\lfloor \sqrt{n-x^2} \rfloor \right) }{6} \right) \\ \end{aligned}

这个式子复杂度已经对了，很多地方可以合并之类的，但是我实在推不动了，之后哪天复活把最简化的结果写在下面（？

代码见下。

#include<stdio.h>
#include<math.h>

typedef long long i64;
typedef double f64;

void assert_val(bool val, char* failed) { if(not val) printf("%s FAILED.\n", failed); }

template <const int mod> struct modint;
template <const int mod> modint<mod> pow(modint<mod> x, int p){
    modint<mod> re(1);
    while(p) {
        if(p&1) re*=x;
        x*=x; p>>=1;
    }
    return re;
}
template <const int mod> struct modint{
    int x;
    modint(int v=0){
        if(v<0 or v>=mod) printf("modint::modint<%d>(%d:int) : Error : %d overflow.\n", mod, v, v);
        this->x = v;
    }
    modint<mod> inv() const { return pow(*this, mod-2); }
};
template <const int mod> bool operator==(modint<mod> a, modint<mod> b) { return a.x==b.x; }
template <const int mod> const modint<mod> operator+(modint<mod> a, modint<mod> b) { return modint<mod>(a.x+b.x>=mod?a.x+b.x-mod:a.x+b.x); }
template <const int mod> const modint<mod> operator-(modint<mod> a) { return a.x==0?0:modint<mod>(mod-a.x); }
template <const int mod> const modint<mod> operator-(modint<mod> a, modint<mod> b) { return a + (-b); }
template <const int mod> const modint<mod> operator*(modint<mod> a, modint<mod> b) { return modint<mod>(1ll * a.x * b.x % mod); }
template <const int mod> const modint<mod> operator/(modint<mod> a, modint<mod> b) { return a * b.inv(); }

template <const int mod> const modint<mod> operator+=(modint<mod>& a, modint<mod> b) { return a = a+b; }
template <const int mod> const modint<mod> operator-=(modint<mod>& a, modint<mod> b) { return a = a-b; }
template <const int mod> const modint<mod> operator*=(modint<mod>& a, modint<mod> b) { return a = a*b; }
template <const int mod> const modint<mod> operator/=(modint<mod>& a, modint<mod> b) { return a = a/b; }

//

const int maxn = 10000007;
const int mod  = 1000000007;

typedef modint<mod> mi;

i64 ni;
mi res;

const mi I1 = mi(1);
const mi I2 = mi(2);
const mi I3 = mi(3);
const mi I4 = mi(4);
const mi I6 = mi(6);
const mi inv2 = mi(2).inv();
const mi inv6 = mi(6).inv();

struct {mi operator[](mi x){return x*(x+mi(1))*(mi(2)*x+mi(1))*inv6;}} pw2sum;

mi pw4sum[maxn], pw6sum[maxn];

mi p2(mi x) { return x*x; }
mi p4(mi x) { return p2(p2(x)); }
mi p6(mi x) { return p2(p2(x)) * p2(x); }

const f64 eps = 0;

// f64 f (f64 x, i64 l) { return x*x - l; }
// f64 df(f64 x) { return 2.*x; }
// i64 sqrt(i64 x) {
    // f64 c = 1;
    // for(int i=1;i<=32;i++)
        // c = c - f(c, x) / df(c);
    // return c + eps;
// }

signed main() {
#ifndef debug_rickyxrc
    freopen("circle.in" , "r", stdin );
    freopen("circle.out", "w", stdout);
#endif

    for(int i=1;i<maxn;i++)
        pw4sum[i] = pw4sum[i-1] + p4(i),
        pw6sum[i] = pw6sum[i-1] + p6(i);

    scanf("%lld", &ni);
    mi n = mi(ni % mod);
    // int qn = sqrt(ni);
    // printf("%d\n", qn);

    // printf("%lld\n", i64(qn + eps) * i64(qn + eps));
    // return 0;

#if 0 // Iter 1
    for(int x=1;x<=sqrt(n);x++)
        for(int y=0;y<=sqrt(n-x*x);y++)
            for(int k=x*x+y*y;k<=n;k++)
                res += mi(k)*mi(n+1-k);
#endif

#if 0 // Iter 2
    for(int x=1;x<=qn;x++)
        for(int y=0;y<=sqrt(ni-x*x);y++)
            for(int ki=x*x+y*y;ki<=ni;ki++) {
                mi k = mi(ki);
                res += (n+I1) * k - k*k;
            }
#endif

#if 0 // Iter 3
    for(int x=1;x<=qn;x++)
        for(int y=0;y<=sqrt(ni-x*x);y++) {
            mi tsum;
            for(int ki=x*x+y*y;ki<=ni;ki++) {
                mi k = mi(ki);
                tsum += k;
            }
            for(int ki=x*x+y*y;ki<=ni;ki++) {
                mi k = mi(ki);
                res -= k*k;
            }
            res += (n+I1)*tsum;
        }
#endif

    // Iter 4
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // mi tsum;
            // tsum = (x*x+y*y+n) * (n-x*x-y*y+I1) * inv2;
            // res += (n+I1)*tsum;
            // for(int ki=1ll*xi*xi+yi*yi;ki<=ni;ki++) {
                // mi k = mi(ki);
                // res -= k*k;
            // }
        // }

    // Iter 5
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // mi tsum;
            // tsum = (x*x+y*y+n) * (n-x*x-y*y+I1) * inv2;
            // res += (n+I1)*tsum;
            // for(int ki=1;ki<1ll*xi*xi+yi*yi;ki++) {
                // mi k = mi(ki);
                // res += k*k;
            // }
            // for(int ki=1;ki<=ni;ki++) {
                // mi k = mi(ki);
                // res -= k*k;
            // }
        // }

    // Iter 6
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += (n+I1) * (x*x+y*y+n) * (n-x*x-y*y+I1) * inv2;
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 7
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += (p2(x)+p2(y)-p4(x)-p4(y)-I2*p2(x)*p2(y)+n*(n+I1));
        // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 8
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * n * (n+I1);
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi y = mi(yi);
            // res += (p2(x)+p2(y)-p4(x)-p4(y)-I2*p2(x)*p2(y));
        // }
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 9
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi y = mi(yi);
            // res += (p2(y)-p4(y));
        // }
        // mi tsum;
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi y = mi(yi);
            // tsum += p2(y);
        // }
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 10
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi y = mi(yi);
            // res += -p4(y);
        // }
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 11
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1], res -= pw2sum[n];
        // }

    // Iter 12
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += pw2sum[x*x+y*y-I1];
        // }
    // }

    // Iter 13
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += (x*x+y*y-I1)*(x*x+y*y)*(I2*x*x+I2*y*y-I1) * inv6;
        // }
    // }

    // Iter 14
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // res += (p4(x)+p4(y)+I2*p2(x)*p2(y)-p2(x)-p2(y))*(I2*x*x+I2*y*y-I1) * inv6;
        // }

    // Iter 15
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
    // mi tsum;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi x = mi(xi), y = mi(yi);
            // tsum += (I2*p6(x) + I2*p6(y)
                // +I6*p2(x)*p4(y) + I6*p4(x)*p2(y)
                // +p2(x)+p2(y)
                // -I6*p2(x)*p2(y)
                // -I3*p4(x)-I3*p4(y));
        // }

    // Iter 16
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        // res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        // mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        // res += -I2*p2(x) * tsum;
    // }
    // res *= (n+I1) * inv2;
    // for(int xi=1;1ll*xi*xi<=ni;xi++)
        // res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
    // mi tsum;
    // for(int xi=1;1ll*xi*xi<=ni;xi++) {
        // mi x = mi(xi);
        // tsum += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (I2*p6(x) + p2(x));
        // for(int yi=0;yi<=sqrt(ni-1ll*xi*xi);yi++) {
            // mi y = mi(yi);
            // tsum += (I2*p6(y)
                // + I6*p2(x)*p4(y) + I6*p4(x)*p2(y)
                // + p2(y)
                // - I6*p2(x)*p2(y)
                // - I3*p4(x)-I3*p4(y));
        // }
    // }

    // Iter 17
    for(int xi=1; 1ll*xi*xi<=ni; xi++) {
        mi x = mi(xi);
        res += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (p2(x) - p4(x) + n * (n+I1));
        res += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        res -= pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        mi tsum = pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        res -= I2 * p2(x) * tsum;
    }
    res *= (n+I1) * inv2;
    for(int xi=1;1ll*xi*xi<=ni;xi++)
        res -= mi(sqrt(ni-1ll*xi*xi) + eps + 1) * pw2sum[n];
    mi tsum;
    for(int xi=1;1ll*xi*xi<=ni;xi++) {
        mi x  = mi(xi);
        tsum += mi(sqrt(ni-1ll*xi*xi) + eps + 1) * (I2*p6(x) + p2(x) - I3*p4(x));
        tsum -= I6 * p2(x) * pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        tsum += I6 * p2(x) * pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        tsum += I2 * pw6sum[int(sqrt(ni-1ll*xi*xi) + eps)];
        tsum += I6 * p4(x) * pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        tsum += pw2sum[sqrt(ni-1ll*xi*xi) + eps];
        tsum -= I3 * pw4sum[int(sqrt(ni-1ll*xi*xi) + eps)];
    }

    //

    res += tsum * inv6;

    res = res * I4 + (n+I1) * (n+I1) * (n) * inv2 - pw2sum[n];

    printf("%d", res.x);
    return 0;
}


在 Rickyxrc's blog 出现的文章，若无特殊注明，均采用 CC BY-NC-SA 4.0 协议共享，也就是转载时需要注明本文章的地址，并且引用本文章的文章也要使用相同的方式共享。