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

“都计算到10^17” ------ 不，我算到10^18。之所以只算到 10^18，那是 unsigned long long 只有64bits，不是算法本身的限制，你换用更高精度的整型就是了。

前两个问题，你自己写一个“大数类”就行了。嫌烦的话，就用别人写好的，比如 GMP (mpz_sqrt)；
最后一个问题，除非你能找到更好的算法，否则计算时间随数据规模增长而增长不是很合理的事吗？我是担心的是，你要 10^99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999……

谢谢各位朋友，已获得程序与数据：
#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

以前我也考虑这个问题，无符号long  long数组。每个元素最大后归0，向上一级元素进一，
思路上好。理解，计算太麻烦。。。。还要考虑这个数组超过了电脑的内存。。。
太麻烦，有没有实用价值呀？？？？