一个计算圆周率任意精度的spigot算法研究

▄︻┻═┳一‥ · 发表于 2011-4-26 10:04:09

这个又让俺想起了，有个叫pade形式有理分式逼近已知函数，非常神奇，对人任意一个函数，用些简单的分式函数逼近到任意精度，比多项式逼近要快多了。

只是俺计算方法这门课学得很不到家，英文也比较烂，有资料却消化不了。

wayne · 发表于 2011-4-26 10:13:55

21# ▄︻┻═┳一‥
是在讨论根号2的近似分数吗

liangbch · 发表于 2011-4-26 10:16:28

用64位整数类型，可将p，q的范围延伸在2^63，倒是double类型的精度只有15到16为10进制数，成为瓶颈。下面的代码给出16位数以内的p和q(用VC6.0编译)

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
typedef unsigned long DWORD;
typedef __int64 INT64;
void print_pq(double f,INT64 limit)
{
INT64 low_p,low_q,high_p,high_q,mid_p,mid_q;
//初始化low_q,low_q,high_p,high_q
low_q=DWORD(floor(f));
high_q=DWORD(ceil(f));
low_p=high_p=1;
mid_p=low_p+high_p;
mid_q=low_q+high_q;
while ( mid_p<limit )
{
printf("%I64u/%I64u=%.16f\n",mid_q,mid_p,double(mid_q)/double(mid_p) );
if ( double(mid_q)/double(mid_p) > f)
{
high_p=mid_p; high_q=mid_q;
}
else
{
low_p=mid_p; low_q=mid_q;
}
mid_p=low_p+high_p;
mid_q=low_q+high_q;
}
}
int main(int argc, char* argv[])
{
print_pq(sqrt(2),10000000000000000I64); //求分子在1000000以内的sqrt(2)有理数逼近
return 0;
}

复制代码

3/2=1.5000000000000000
4/3=1.3333333333333333
7/5=1.3999999999999999
10/7=1.4285714285714286
17/12=1.4166666666666667
24/17=1.4117647058823530
41/29=1.4137931034482758
58/41=1.4146341463414633
99/70=1.4142857142857144
140/99=1.4141414141414141
239/169=1.4142011834319526
338/239=1.4142259414225942
577/408=1.4142156862745099
816/577=1.4142114384748701
1393/985=1.4142131979695431
1970/1393=1.4142139267767408
3363/2378=1.4142136248948696
4756/3363=1.4142134998513232
8119/5741=1.4142135516460548
11482/8119=1.4142135731001355
19601/13860=1.4142135642135643
27720/19601=1.4142135605326258
47321/33461=1.4142135620573204
66922/47321=1.4142135626888697
114243/80782=1.4142135624272734
161564/114243=1.4142135623189167
275807/195025=1.4142135623637995
390050/275807=1.4142135623823906
665857/470832=1.4142135623746899
941664/665857=1.4142135623715002
1607521/1136689=1.4142135623728214
2273378/1607521=1.4142135623733687
3880899/2744210=1.4142135623731420
5488420/3880899=1.4142135623730481
9369319/6625109=1.4142135623730869
13250218/9369319=1.4142135623731031
22619537/15994428=1.4142135623730965
31988856/22619537=1.4142135623730936
54608393/38613965=1.4142135623730947
77227930/54608393=1.4142135623730954
131836323/93222358=1.4142135623730951
209064253/147830751=1.4142135623730951
286292183/202439144=1.4142135623730951
363520113/257047537=1.4142135623730951
440748043/311655930=1.4142135623730951
517975973/366264323=1.4142135623730951
595203903/420872716=1.4142135623730951
672431833/475481109=1.4142135623730951
749659763/530089502=1.4142135623730951
826887693/584697895=1.4142135623730951
904115623/639306288=1.4142135623730954
1731003316/1224004183=1.4142135623730951
2635118939/1863310471=1.4142135623730951
3539234562/2502616759=1.4142135623730954
6174353501/4365927230=1.4142135623730954
8809472440/6229237701=1.4142135623730954
11444591379/8092548172=1.4142135623730954
14079710318/9955858643=1.4142135623730951
25524301697/18048406815=1.4142135623730951
36968893076/26140954987=1.4142135623730954
62493194773/44189361802=1.4142135623730954
88017496470/62237768617=1.4142135623730954
113541798167/80286175432=1.4142135623730954
139066099864/98334582247=1.4142135623730954
164590401561/116382989062=1.4142135623730954
190114703258/134431395877=1.4142135623730954
215639004955/152479802692=1.4142135623730954
241163306652/170528209507=1.4142135623730954
266687608349/188576616322=1.4142135623730954
292211910046/206625023137=1.4142135623730954
317736211743/224673429952=1.4142135623730954
343260513440/242721836767=1.4142135623730954
368784815137/260770243582=1.4142135623730954
394309116834/278818650397=1.4142135623730954
419833418531/296867057212=1.4142135623730954
445357720228/314915464027=1.4142135623730951
865191138759/611782521239=1.4142135623730954
1310548858987/926697985266=1.4142135623730951
2175739997746/1538480506505=1.4142135623730951
3040931136505/2150263027744=1.4142135623730951
3906122275264/2762045548983=1.4142135623730951
4771313414023/3373828070222=1.4142135623730951
5636504552782/3985610591461=1.4142135623730951
6501695691541/4597393112700=1.4142135623730954
12138200244323/8583003704161=1.4142135623730954
17774704797105/12568614295622=1.4142135623730954
23411209349887/16554224887083=1.4142135623730954
29047713902669/20539835478544=1.4142135623730954
34684218455451/24525446070005=1.4142135623730951
63731932358120/45065281548549=1.4142135623730951
92779646260789/65605117027093=1.4142135623730951
121827360163458/86144952505637=1.4142135623730951
150875074066127/106684787984181=1.4142135623730951
179922787968796/127224623462725=1.4142135623730954
330797862034923/233909411446906=1.4142135623730954
481672936101050/340594199431087=1.4142135623730954
632548010167177/447278987415268=1.4142135623730951
1114220946268227/787873186846355=1.4142135623730951
1595893882369277/1128467386277442=1.4142135623730954
2710114828637504/1916340573123797=1.4142135623730954
3824335774905731/2704213759970152=1.4142135623730954
4938556721173958/3492086946816507=1.4142135623730954
6052777667442185/4279960133662862=1.4142135623730954
7166998613710412/5067833320509217=1.4142135623730954
8281219559978639/5855706507355572=1.4142135623730954
9395440506246866/6643579694201927=1.4142135623730954
10509661452515093/7431452881048282=1.4142135623730954
11623882398783320/8219326067894637=1.4142135623730954
12738103345051547/9007199254740992=1.4142135623730954
13852324291319774/9795072441587347=1.4142135623730951

复制代码

liangbch · 发表于 2011-4-26 10:28:56

本帖最后由 liangbch 于 2011-4-26 10:31 编辑

22# wayne

16# 的目地是求N,a,p,q的值。
我的代码也是求p,q 的值，只不过没有给出N,a的值而已. ▄︻┻═┳一‥的代码大概是穷举法，速度很慢，当N很大时，慢的无法接受。我的是迭代法，速度很快。

by the way, 我的代码适用于m=2的情况

▄︻┻═┳一‥ · 发表于 2011-4-26 10:29:36

21# ▄︻┻═┳一‥
是在讨论根号2的近似分数吗
wayne 发表于 2011-4-26 10:13

不是的，版主说用分数q/p逼近一个浮点数，让我想到了连分数的最佳有理逼近特性。。

也不知liangbch版主为什么不用连分数逼近？？

liangbch · 发表于 2011-4-26 10:33:30

25# ▄︻┻═┳一‥

从本质上，18楼表中的p和q, p/q就是 sqrt(2)的一个分数逼近。

▄︻┻═┳一‥ · 发表于 2011-4-26 10:39:57

大家好像对那个不定方程里的p,q有误解

它只是那个像秦九绍多项式的系数，和$root{m}{N}$没什么联系的

看#9楼的那个级数就知道了。

▄︻┻═┳一‥ · 发表于 2011-4-26 10:50:06

如果仅仅是逼近一个二次的无理数，那么连分数里有相关定理：

拉格朗日定理：任何一个二次无理数都有一个无限连分数展式，这个展示从点往后是循环的。
进一步：sqrt(N)的展示循环节总是形如[a1, (a2, a3, ... ,an, 2*a1)]，这里（a2, a3, ... ,an, 2*a1）是展示的循环节；它的最末位就是a1的二倍；

▄︻┻═┳一‥ · 发表于 2011-4-26 11:22:48

当m=2时
|b|<10,100以内整数解：

#-----+-----+-----+-----+-----+
# N | a | b | p | q |
#-----+-----+-----+-----+-----+
# 2| 8| -1| 2| 3|
# 2| 9| 1| 3| 4|
# 2| 9| 7| 3| 2|
# 2| 18| -7| 3| 5|
# 2| 25| 7| 5| 6|
# 2| 25| -7| 5| 8|
# 2| 32| 7| 4| 5|
# 2| 49| -1| 7| 10|
# 2| 50| 1| 5| 7|
# 3| 4| 1| 2| 3|
# 3| 25| -2| 5| 9|
# 3| 27| 2| 3| 5|
# 3| 48| -1| 4| 7|
# 3| 49| 1| 7| 12|
# 5| 4| -1| 2| 5|
# 5| 9| 4| 3| 5|
# 5| 45| -4| 3| 7|
# 5| 49| 4| 7| 15|
# 5| 80| -1| 4| 9|
# 5| 81| 1| 9| 20|
# 6| 8| 5| 2| 3|
# 6| 24| -1| 2| 5|
# 6| 25| 1| 5| 12|
# 6| 27| -5| 3| 8|
# 6| 32| 5| 4| 9|
# 6| 49| -5| 7| 18|
# 6| 54| 5| 3| 7|
# 7| 4| -3| 2| 7|
# 7| 9| 2| 3| 7|
# 7| 16| 9| 4| 7|
# 7| 25| -3| 5| 14|
# 7| 28| 3| 2| 5|
# 7| 63| -1| 3| 8|
# 7| 64| 1| 8| 21|
# 8| 9| 1| 3| 8|
# 8| 9| 7| 3| 4|
# 8| 18| -7| 3| 10|
# 8| 25| 7| 5| 12|
# 8| 25| -7| 5| 16|
# 8| 32| 7| 2| 5|
# 8| 49| -1| 7| 20|
# 8| 50| 1| 5| 14|
# 10| 8| 3| 2| 5|
# 10| 9| -1| 3| 10|
# 10| 40| -9| 2| 7|
# 10| 49| 9| 7| 20|
# 11| 4| -7| 2| 11|
# 11| 9| -2| 3| 11|
# 11| 16| 5| 4| 11|
# 11| 44| -5| 2| 7|
# 11| 49| 5| 7| 22|
# 11| 99| -1| 3| 10|
# 12| 25| -2| 5| 18|
# 12| 27| 2| 3| 10|
# 12| 48| -1| 2| 7|
# 12| 49| 1| 7| 24|
# 13| 4| -9| 2| 13|
# 13| 9| -4| 3| 13|
# 13| 16| 3| 4| 13|
# 13| 49| -3| 7| 26|
# 13| 52| 3| 2| 7|
# 14| 8| 1| 2| 7|
# 14| 9| -5| 3| 14|
# 15| 12| 7| 2| 5|
# 15| 16| 1| 4| 15|
# 15| 20| -7| 2| 9|
# 15| 27| 7| 3| 10|
# 17| 9| -8| 3| 17|
# 17| 16| -1| 4| 17|
# 17| 25| 8| 5| 17|
# 18| 8| -1| 2| 9|
# 18| 8| 7| 2| 3|
# 18| 25| 7| 5| 18|
# 18| 25| -7| 5| 24|
# 18| 32| 7| 4| 15|
# 18| 49| -1| 7| 30|
# 18| 50| 1| 5| 21|
# 19| 16| -3| 4| 19|
# 19| 25| 6| 5| 19|
# 19| 76| -5| 2| 9|
# 19| 81| 5| 9| 38|
# 20| 9| 4| 3| 10|
# 20| 45| -4| 3| 14|
# 20| 49| 4| 7| 30|
# 20| 80| -1| 2| 9|
# 20| 81| 1| 9| 40|
# 21| 12| 5| 2| 7|
# 21| 16| -5| 4| 21|
# 21| 25| 4| 5| 21|
# 21| 27| -1| 3| 14|
# 21| 28| 1| 2| 9|
# 22| 8| -3| 2| 11|
# 22| 18| 7| 3| 11|
# 22| 25| 3| 5| 22|
# 22| 81| -7| 9| 44|
# 22| 88| 7| 2| 9|
# 22| 98| -1| 7| 33|
# 22| 99| 1| 3| 14|
# 23| 16| -7| 4| 23|
# 23| 25| 2| 5| 23|
# 24| 25| 1| 5| 24|
# 24| 27| -5| 3| 16|
# 24| 32| 5| 2| 9|
# 24| 49| -5| 7| 36|
# 24| 54| 5| 3| 14|
# 26| 8| -5| 2| 13|
# 26| 18| 5| 3| 13|
# 26| 25| -1| 5| 26|
# 27| 4| 1| 2| 9|
# 27| 25| -2| 5| 27|
# 27| 48| -1| 4| 21|
# 27| 49| 1| 7| 36|
# 28| 9| 2| 3| 14|
# 28| 16| 9| 2| 7|
# 28| 25| -3| 5| 28|
# 28| 63| -1| 3| 16|
# 28| 64| 1| 4| 21|
# 29| 25| -4| 5| 29|
# 29| 36| 7| 6| 29|
# 30| 8| -7| 2| 15|
# 31| 25| -6| 5| 31|
# 31| 36| 5| 6| 31|
# 32| 9| 1| 3| 16|
# 32| 9| 7| 3| 8|
# 32| 18| -7| 3| 20|
# 32| 25| 7| 5| 24|
# 32| 25| -7| 5| 32|
# 32| 49| -1| 7| 40|
# 32| 50| 1| 5| 28|
# 33| 12| 1| 2| 11|
# 33| 25| -8| 5| 33|
# 34| 8| -9| 2| 17|
# 34| 18| 1| 3| 17|
# 34| 25| -9| 5| 34|
# 35| 36| 1| 6| 35|
# 37| 36| -1| 6| 37|
# 38| 18| -1| 3| 19|
# 39| 12| -1| 2| 13|
# 40| 9| -1| 3| 20|
# 40| 49| 9| 7| 40|
# 41| 36| -5| 6| 41|
# 41| 49| 8| 7| 41|
# 43| 36| -7| 6| 43|
# 43| 49| 6| 7| 43|
# 44| 9| -2| 3| 22|
# 44| 16| 5| 2| 11|
# 44| 49| 5| 7| 44|
# 44| 99| -1| 3| 20|
# 45| 4| -1| 2| 15|
# 45| 49| 4| 7| 45|
# 45| 80| -1| 4| 27|
# 45| 81| 1| 3| 20|
# 46| 18| -5| 3| 23|
# 46| 32| 9| 4| 23|
# 46| 49| 3| 7| 46|
# 47| 49| 2| 7| 47|
# 48| 25| -2| 5| 36|
# 48| 27| 2| 3| 20|
# 48| 49| 1| 7| 48|
# 50| 8| -1| 2| 15|
# 50| 8| 7| 2| 5|
# 50| 9| 1| 3| 20|
# 50| 9| 7| 3| 10|
# 50| 18| -7| 3| 25|
# 50| 32| 7| 4| 25|
# 50| 49| -1| 7| 50|
# 51| 12| -5| 2| 17|
# 51| 49| -2| 7| 51|
# 51| 68| -7| 2| 15|
# 51| 75| 7| 5| 34|
# 52| 9| -4| 3| 26|
# 52| 16| 3| 2| 13|
# 52| 49| -3| 7| 52|
# 53| 49| -4| 7| 53|
# 54| 8| 5| 2| 9|
# 54| 24| -1| 2| 15|
# 54| 25| 1| 5| 36|
# 54| 32| 5| 4| 27|
# 54| 49| -5| 7| 54|
# 55| 20| 9| 2| 11|
# 55| 44| -1| 2| 15|
# 55| 45| 1| 3| 22|
# 55| 49| -6| 7| 55|
# 55| 64| 9| 8| 55|
# 56| 9| -5| 3| 28|
# 57| 12| -7| 2| 19|
# 57| 27| 8| 3| 19|
# 57| 49| -8| 7| 57|
# 57| 64| 7| 8| 57|
# 57| 75| -1| 5| 38|
# 57| 76| 1| 2| 15|
# 58| 32| 3| 4| 29|
# 58| 49| -9| 7| 58|
# 59| 64| 5| 8| 59|
# 60| 16| 1| 2| 15|
# 60| 27| 7| 3| 20|
# 61| 64| 3| 8| 61|
# 62| 32| 1| 4| 31|
# 63| 4| -3| 2| 21|
# 63| 16| 9| 4| 21|
# 63| 25| -3| 5| 42|
# 63| 28| 3| 2| 15|
# 63| 64| 1| 8| 63|
# 65| 20| 7| 2| 13|
# 65| 45| -7| 3| 26|
# 65| 52| 7| 2| 15|
# 65| 64| -1| 8| 65|
# 66| 27| 5| 3| 22|
# 66| 32| -1| 4| 33|
# 67| 64| -3| 8| 67|
# 68| 9| -8| 3| 34|
# 68| 16| -1| 2| 17|
# 68| 25| 8| 5| 34|
# 69| 27| 4| 3| 23|
# 69| 64| -5| 8| 69|
# 70| 32| -3| 4| 35|
# 71| 64| -7| 8| 71|
# 72| 25| 7| 5| 36|
# 72| 25| -7| 5| 48|
# 72| 32| 7| 2| 15|
# 72| 49| -1| 7| 60|
# 72| 50| 1| 5| 42|
# 73| 64| -9| 8| 73|
# 73| 81| 8| 9| 73|
# 74| 32| -5| 4| 37|
# 74| 81| 7| 9| 74|
# 75| 4| 1| 2| 15|
# 75| 27| 2| 3| 25|
# 75| 48| -1| 4| 35|
# 75| 49| 1| 7| 60|
# 76| 16| -3| 2| 19|
# 76| 25| 6| 5| 38|
# 76| 81| 5| 9| 76|
# 77| 81| 4| 9| 77|
# 78| 27| 1| 3| 26|
# 78| 32| -7| 4| 39|
# 79| 81| 2| 9| 79|
# 80| 9| 4| 3| 20|
# 80| 45| -4| 3| 28|
# 80| 49| 4| 7| 60|
# 80| 81| 1| 9| 80|
# 82| 32| -9| 4| 41|
# 82| 50| 9| 5| 41|
# 82| 81| -1| 9| 82|
# 83| 81| -2| 9| 83|
# 84| 16| -5| 2| 21|
# 84| 25| 4| 5| 42|
# 84| 27| -1| 3| 28|
# 85| 20| 3| 2| 17|
# 85| 81| -4| 9| 85|
# 86| 50| 7| 5| 43|
# 86| 81| -5| 9| 86|
# 87| 27| -2| 3| 29|
# 88| 18| 7| 3| 22|
# 88| 25| 3| 5| 44|
# 88| 81| -7| 9| 88|
# 88| 98| -1| 7| 66|
# 88| 99| 1| 3| 28|
# 89| 81| -8| 9| 89|
# 90| 8| 3| 2| 15|
# 90| 40| -9| 2| 21|
# 90| 49| 9| 7| 60|
# 92| 16| -7| 2| 23|
# 92| 25| 2| 5| 46|
# 93| 27| -4| 3| 31|
# 94| 50| 3| 5| 47|
# 95| 20| 1| 2| 19|
# 96| 25| 1| 5| 48|
# 96| 27| -5| 3| 32|
# 96| 49| -5| 7| 72|
# 96| 54| 5| 3| 28|
# 98| 8| -1| 2| 21|
# 98| 8| 7| 2| 7|
# 98| 9| 1| 3| 28|
# 98| 9| 7| 3| 14|
# 98| 18| -7| 3| 35|
# 98| 25| 7| 5| 42|
# 98| 25| -7| 5| 56|
# 98| 32| 7| 4| 35|
# 98| 50| 1| 5| 49|
# 99| 4| -7| 2| 33|
# 99| 16| 5| 4| 33|
# 99| 44| -5| 2| 21|
# 99| 49| 5| 7| 66|

复制代码

▄︻┻═┳一‥ · 发表于 2011-4-26 11:24:39

看来，如果放宽b的限制，让b可以为负数，这个级数适用范围还是蛮宽的，，，，

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 一个计算圆周率任意精度的spigot算法研究