输出园半径R的园内整点个数，园周率精确到40位有效值，

版主好！
半径R的园内 u^2 + v^2 ≤x= R^2，整点个数f(x)，已知数据如下：
R        f(x)
10^1     37
10^2     317
10^3     3149
10^4     31417
10^5     314197
10^6     3141549
10^7     31416025
10^8     314159053
10^9     3141592409
10^10    31415925457
10^11    314159264013
10^12    3141592649625
10^13    31415926532017
10^14    314159265350589
10^15    3141592653588533
10^16    31415926535867961
10^17    314159265358987341
希望计算更大x值的f(x)值，

[此贴子已经被作者于2022-5-21 23:40编辑过]

定义：间距是1的水平平行线簇，与间距是1的垂直平行线簇的交点，称为整点；或：u,v都是整数构成的二维坐标点（u,v），称为整点。
可以推广到多维情况；
设本高斯园的园心位于坐标原点（0，0），只需要计算并累加园内 2R 排（或2R列）的整点个数，某u值的水平整点个数f(x,u)=2（x-v^2）^0.5,

f(x)=Σf(x,u)=2（x-v^2）^0.5
其中，u值的求和范围（0，[x^0.5]）

谢谢wp231957指正：统计并累加园内整点个数确定与园周率无关，理论值πx，才与园周率有关

[此贴子已经被作者于2022-5-22 08:45编辑过]

谢谢rjsp编程计算！！！与10多年前得到的数据（没程序）一样，都计算到10^17，可能两种方法相同，并遇到相同的问题，不能计算更大数据，素数定理计算到10^28，不知道用什么方法：
（https://bbs.emath.）
rjsp先生提出的三个问题，能解决嘛？

如果半径为10：则：x=100,理论个数=πx≈314，上面编程计算的实际个数=317，显然，x越大，实际个数越接近理论个数πx

理论上，数学期望是πx ，最大误差约是2.31205*[(x)^0.25]*(lnlnx)^0.5，在10^17内，实际数据是符合理论分析解

谢谢各位朋友，已获得程序与数据：
#include <cstdio>
#include <cmath>
#include <cstdint>

#include <utility>
#include <stack>
#include <functional>

using namespace std;
using i64 = int64_t;

i64 solve_fast(i64 N) {
auto inside = [N] (i64 x, i64 y) {
return x * x + y * y <= N;
};
auto cut = [] (i64 x, i64 y, int dx1, int dy1) {
return dx1 * x >= dy1 * y;
};

const i64 v = sqrtl(N / 2), w = sqrtl(N);
i64 x = v;
i64 y = i64(sqrtl(max<i64>(0, N - (v + 1) * (v + 1)))) + 1;

auto stac = stack< pair<int, int> >({{0, 1}, {1, 1}});

i64 ret = 0;
while (1) {
int dx1, dy1; tie(dx1, dy1) = stac.top(); stac.pop();
while (inside(x + dx1, y - dy1)) {
x += dx1; y -= dy1;

ret += i64(dx1) * (y - 1)
+ ((i64(dx1 + 1) * (dy1 + 1)) >> 1) - dy1;
}

int dx2 = dx1, dy2 = dy1;
while (!stac.empty()) {
tie(dx1, dy1) = stac.top();
if (inside(x + dx1, y - dy1)) break;
stac.pop();
dx2 = dx1, dy2 = dy1;
}
if (stac.empty()) break;

while (1) {
int dx12 = dx1 + dx2, dy12 = dy1 + dy2;
if (inside(x + dx12, y - dy12)) {
stac.emplace(dx1 = dx12, dy1 = dy12);
} else {
if (cut(x + dx12, y - dy12, dx1, dy1)) break;
dx2 = dx12, dy2 = dy12;
}
}
}
ret = ret * 2 + i64(v) * v;
ret = ret * 4 + 4 * i64(w) + 1;
return ret;
}

i64 solve_naive(i64 N) {
int v = (int) sqrtl(N);
i64 ret = 0;
for (int i = 1; i <= v; ++i) {
ret += (int) sqrtl(N - i64(i) * i);
}
return 4 * (ret + v) + 1;
}

int main() {
i64 N = 1e18;
//printf("%llu\n", solve_naive(N));
printf("%llu\n", solve_fast(N));
return 0;
}



这个cpp代码可以算到10^18，再大范围就需要修改数据类型了
0 5
1 37
2 317
3 3149
4 31417
5 314197
6 3141549
7 31416025
8 314159053
9 3141592409
10 31415925457
11 314159264013
12 3141592649625
13 31415926532017
14 314159265350589
15 3141592653588533
16 31415926535867961
17 314159265358987341
18 3141592653589764829
19 31415926535897744669
20 314159265358978759661
21 3141592653589792630933
22 31415926535897931085161
23 314159265358979322639853
24 3141592653589793234680617
25 31415926535897932384615349
26 314159265358979323823745421
27 3141592653589793238428435569
28 31415926535897932384568540625
29 314159265358979323846212602093
30 3141592653589793238462579472373
31 31415926535897932384626459376945
32 314159265358979323846263865968245
33 3141592653589793238462643289640533
34 31415926535897932384626432234171745
35 314159265358979323846264338399627025
36 3141592653589793238462643379445627833