Link
考虑点((x,y))的贡献，设(k=x^2+y^2)，那么贡献为：
[egin{aligned}
&sumlimits_{i=k}^nsumlimits_{j=k}^ij\
=&sumlimits_{i=k}^n{i+1choose 2}-{kchoose 2}\
=&{n+2choose 3}-{kchoose 3}-(n-k-1){kchoose 2}
end{aligned}
]
枚举(x)，首先可以得到(y)的取值范围。
同时我们知道(k)是一个关于(y)的二次多项式，因此贡献是关于(y)的六次多项式，再套上一个(sum)就是七次多项式。
预处理幂和即可(O(1))计算，总时间复杂度为(O(sqrt n))。
#include<cmath>
#include<cstdio>
using i64=long long;
const int N=1000007,P=1000000007;
i64 n,m,md,a[N],b[N],c[N],ans;
int main()
{
scanf("%lld",&n),m=ceil(sqrt(n)+1e-8)-1,md=n%P;
for(i64 i=1,t;i<=m;++i) t=i*i%P,a[i]=(a[i-1]+t)%P,b[i]=(b[i-1]+t*t)%P,c[i]=(c[i-1]+t*t%P*t)%P;
for(i64 i=1,j=m,k;i<=m;++i)
{
for(;i*i+j*j>n;--j);
k=i*i%P,ans=(ans+(1+md-k)*(2+md-k)%P*(md+k+k)%P*(j+1)+(4+3*md-12*k%P-6*md*k%P+6*k*k+P+P)%P*a[j]+(6*k-3*md%P-6+P+P)%P*b[j]+2*c[j])%P;
}
ans=2*((P+1)/3)*ans%P,ans=(ans+md*(md+1)%P*(md+2)%P*((P+1)/6))%P,printf("%lld",ans);
}