有偿求助：编写统计园内整点个数的程序 - 编程论坛

#2

rjsp2022-05-20 08:32

有谁能听懂他在说什么吗？
是不是说求满足 x^2 + y^2 <= r^2 条件的所有(x,y)整数坐标的个数？如果这样的话，跟“圆周率”有什么屁关系？
还有“精确到40位有效值”，float才7.22位，double才15.95位，long double才19.26位。如果要40位精确度，那么这种浮点数的有效尾数就要133位，连非标准的float128都无法满足你的要求。

#3

njzzyy2022-05-21 22:33

版主好！
半径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编辑过]

#4

wp2319572022-05-22 06:00

回复 3楼 njzzyy

这是求π啊，你还是没说啥是整点啊，管说有多少个谁知道你那数是咋来的啊
再说了10的17次幂已经超范围了，普通计算机算不了太大的数

#5

njzzyy2022-05-22 08:39

定义：间距是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编辑过]

#6

rjsp2022-05-22 13:54

我觉得有三个难点
a. 用什么类型来存储更大的数（可以写一个大数类）
b. 怎么快速开根号（可以用华罗庚快速开根号法）
c. 计算量是0.5倍指数上涨的（除非有更好的算法，但几乎不可能了）

程序代码：

#include <math.h>
#include <assert.h>

// 对一个整数开根号，这是效率的关键所在，但我使用了最烂的方法，只为了验证算法的正确性
unsigned long long sqrtint( unsigned long long n )
{
    unsigned long long r = (unsigned long long)sqrt( n*1.0 );
    // 因为sqrt精度不足，所以需要校正
    if( r*r > n )
        for( --r; r*r>n; --r );
    else if( r*r < n )
        for( ; (r+1)*(r+1)<=n; ++r );
    return r;
}

// 求 x^2 + y^2 <= RR 整数解的数目, 结果大约是 PI*RR
unsigned long long foo( unsigned long long RR )
{
    unsigned long long r = sqrtint( RR );
    unsigned long long w = sqrtint( RR/2 );

    unsigned long long t = 0;
    for( unsigned long long i=w+1; i<=r; ++i )
        t += sqrtint( RR - i*i );

    unsigned long long result = (2*w+1)*(2*w+1) + 4*(r-w) + 8*t;
    // 上述表达式可能产生数据溢出，另一种方法是在函数入口就限定 RR*PI必须小于unsigned long long所能表达的最大数
    assert( (2*w+1)*(2*w+1)/(2*w+1) == (2*w+1) );
    assert( result > 4*(r-w) + 8*t );
    return result;
}

#include <stdio.h>

int main( void )
{
    assert( foo(10) == 37 );
    assert( foo(100) == 317 );
    assert( foo(1000) == 3149 );
    assert( foo(10000) == 31417 );
    assert( foo(100000) == 314197 );
    assert( foo(1000000) == 3141549 );
    assert( foo(10000000) == 31416025 );
    assert( foo(100000000) == 314159053 );
    assert( foo(1000000000) == 3141592409 );
    assert( foo(10000000000) == 31415925457 );
    assert( foo(100000000000) == 314159264013 );
    assert( foo(1000000000000) == 3141592649625 );
    assert( foo(10000000000000) == 31415926532017 );
    assert( foo(100000000000000) == 314159265350589 );
    assert( foo(1000000000000000) == 3141592653588533 );
    assert( foo(10000000000000000) == 31415926535867961 );
    assert( foo(100000000000000000) == 314159265358987341 );
    printf( "f(10^18) = %llu\n", foo(1000000000000000000) );
    printf( "f(10^19) = %llu\n", foo(10000000000000000000) );
}

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

#7

njzzyy2022-05-22 21:13

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

#8

jklqwe1112022-05-22 23:28

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

无法理解这些数据，如果半径为10，那么，互相垂直的两条直径上就有41个整点，全部计算出来，差距很大。也许我理解有误，楼主能否解惑。
计算公式大概应该是这样 ∑2*√(r*r-v*v)+1 v ∈[ -r,r] 公式正确性有待证明。

#9

jklqwe1112022-05-23 01:08

计算公式 r为奇数: f(r)=((r/2)*(r+1)+(r+1)/2)*4+1
r为偶数: f(r)=((r/2)*(r+1))*4+1

f(r)=4*(1+...+r)+1

#10

njzzyy2022-05-23 08:40

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

#11

njzzyy2022-05-23 08:46

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

#12

rjsp2022-05-23 10:17

以下是引用njzzyy在2022-5-22 21:13:04的发言：

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

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

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

#13

njzzyy2022-05-23 15:31

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

#14

追梦人zmrghy2022-05-25 12:55

回复 7楼 njzzyy

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