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

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

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

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

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

我觉得有三个难点
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编辑过]

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

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] 公式正确性有待证明。

计算公式 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：则：x=100,理论个数=πx≈314，上面编程计算的实际个数=317，显然，x越大，实际个数越接近理论个数πx