▄︻┻═┳一‥ 发表于 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 http://bbs.emath.ac.cn/images/common/back.gif

不是的,版主说用分数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)的展示循环节总是形如,这里(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可以为负数,这个级数适用范围还是蛮宽的,,,,
页: 1 2 [3] 4 5 6 7 8
查看完整版本: 一个计算圆周率任意精度的spigot算法研究