算法挑战：为浮子生成持续的分数
https://stackoverflow.com/questions/4637967
-
08-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian题
(编辑: ：回答脾气暴躁的评论，不，这不是作业。我正在研究音高检测，采用一系列潜在的谐波峰，并试图为基本频率构建候选者。因此，这实际上是一个非常实用的问题。）
考虑（例如）PI的最佳分数近似值，通过增加分母排序：3/1、22/7、355/113，...
挑战：创建一个 整齐的 C算法将生成给定浮子的n'th商近似A/B，也返回差异。
ColcBestFrac（float frac，int n，int * a，int * b，float * err）{...}
我相信最好的技术是 继续分数
拿走PI的分数部分，您会得到3
现在，其余为0.14159 ... = 1/7.06251 ..
因此，下一个最好的理性是3 + 1/7 = 22/7
从7.06251拿走7，您将获得0.06251 ..大约1/15.99659 ..
称其为16，然后下一个最佳近似是
3 + 1/(7 + 1/16) = 355/113
但是，这远非微不足道的转换为干净的C代码。如果我得到一些整理的话，我会发布。同时，某人可能会喜欢它作为Brainteaser。
解决方案
由于您要求这个作为答案而不是评论。
对于任何实际数字，其持续分数的收敛p [k]/q [k]始终是最好的理性近似值，但不是 全部 最好的理性近似值。要获取所有这些 (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) 对于一些整数n≥1。取n = a [k+2]给出p [k+2]/q [k+2]，而整数n是从地板（a [k+2]/2）或天花板（a [a [a [ k+2]/2），到[k+2]。这也提到 在维基百科.
近似π
π的持续分数为[3; 7、15、1、292、1、1、1、2、1、3、1、14、2…]（序列 A001203在OEIS），收敛顺序为3/1、22/7、333/106、355/113、103993/33102…（（A002485/A002486），最佳近似顺序为3/1、13/4、16/5、19/6、22/7、179/57…（A063674/A063673).
因此，该算法说π= [3; 7、15、1、292、1、1，…
3/1 = [3]
13/4 = [3; 4]
16/5 = [3; 5]
19/6 = [3; 6]
22/7 = [3; 7]
179/57 = [3; 7, 8]
201/64 = [3; 7, 9]
223/71 = [3; 7, 10]
245/78 = [3; 7, 11]
267/85 = [3; 7, 12]
289/92 = [3; 7, 13]
311/99 = [3; 7, 14]
333/106 = [3; 7, 15]
355/113 = [3; 7, 15, 1]
52163/16604 = [3; 7, 15, 1, 146]
52518/16717 = [3; 7, 15, 1, 147]
… (all the fractions from [3; 7, 15, 1, 148] to [3; 7, 15, 1, 291])…
103993/33102 = [3; 7, 15, 1, 292]
104348/33215 = [3; 7, 15, 1, 292, 1]
...
程序
这是一个C程序，该程序给出了一个正实数，生成其持续分数，收敛和最佳理性近似序列。功能 find_cf 找到持续的分数（将术语放在[]和Q []中的汇编中，请原谅全局变量），函数 all_best 打印所有最佳理性近似值。
#include <math.h>
#include <stdio.h>
#include <assert.h>
// number of terms in continued fraction.
// 15 is the max without precision errors for M_PI
#define MAX 15
#define eps 1e-9
long p[MAX], q[MAX], a[MAX], len;
void find_cf(double x) {
int i;
//The first two convergents are 0/1 and 1/0
p[0] = 0; q[0] = 1;
p[1] = 1; q[1] = 0;
//The rest of the convergents (and continued fraction)
for(i=2; i<MAX; ++i) {
a[i] = lrint(floor(x));
p[i] = a[i]*p[i-1] + p[i-2];
q[i] = a[i]*q[i-1] + q[i-2];
printf("%ld: %ld/%ld\n", a[i], p[i], q[i]);
len = i;
if(fabs(x-a[i])<eps) return;
x = 1.0/(x - a[i]);
}
}
void all_best(double x) {
find_cf(x); printf("\n");
int i, n; long cp, cq;
for(i=2; i<len; ++i) {
//Test n = a[i+1]/2. Enough to test only when a[i+1] is even, actually...
n = a[i+1]/2; cp = n*p[i]+p[i-1]; cq = n*q[i]+q[i-1];
if(fabs(x-(double)cp/cq) < fabs(x-(double)p[i]/q[i]))
printf("%ld/%ld, ", cp, cq);
//And print all the rest, no need to test
for(n = (a[i+1]+2)/2; n<=a[i+1]; ++n) {
printf("%ld/%ld, ", n*p[i]+p[i-1], n*q[i]+q[i-1]);
}
}
}
int main(int argc, char **argv) {
double x;
if(argc==1) { x = M_PI; } else { sscanf(argv[1], "%lf", &x); }
assert(x>0); printf("%.15lf\n\n", x);
all_best(x); printf("\n");
return 0;
}
例子
对于π，这是该程序的输出，在大约0.003秒内（即，它比循环浏览所有可能的分母！），可读性：可读性：
% ./a.out
3.141592653589793
3: 3/1
7: 22/7
15: 333/106
1: 355/113
292: 103993/33102
1: 104348/33215
1: 208341/66317
1: 312689/99532
2: 833719/265381
1: 1146408/364913
3: 4272943/1360120
1: 5419351/1725033
14: 80143857/25510582
13/4, 16/5, 19/6, 22/7, 179/57, 201/64, 223/71, 245/78, 267/85, 289/92, 311/99,
333/106, 355/113, 52163/16604, 52518/16717, 52873/16830, 53228/16943, 53583/17056,
53938/17169, 54293/17282, 54648/17395, 55003/17508, 55358/17621, 55713/17734,
56068/17847, 56423/17960, 56778/18073, 57133/18186, 57488/18299, 57843/18412,
58198/18525, 58553/18638, 58908/18751, 59263/18864, 59618/18977, 59973/19090,
60328/19203, 60683/19316, 61038/19429, 61393/19542, 61748/19655, 62103/19768,
62458/19881, 62813/19994, 63168/20107, 63523/20220, 63878/20333, 64233/20446,
64588/20559, 64943/20672, 65298/20785, 65653/20898, 66008/21011, 66363/21124,
66718/21237, 67073/21350, 67428/21463, 67783/21576, 68138/21689, 68493/21802,
68848/21915, 69203/22028, 69558/22141, 69913/22254, 70268/22367, 70623/22480,
70978/22593, 71333/22706, 71688/22819, 72043/22932, 72398/23045, 72753/23158,
73108/23271, 73463/23384, 73818/23497, 74173/23610, 74528/23723, 74883/23836,
75238/23949, 75593/24062, 75948/24175, 76303/24288, 76658/24401, 77013/24514,
77368/24627, 77723/24740, 78078/24853, 78433/24966, 78788/25079, 79143/25192,
79498/25305, 79853/25418, 80208/25531, 80563/25644, 80918/25757, 81273/25870,
81628/25983, 81983/26096, 82338/26209, 82693/26322, 83048/26435, 83403/26548,
83758/26661, 84113/26774, 84468/26887, 84823/27000, 85178/27113, 85533/27226,
85888/27339, 86243/27452, 86598/27565, 86953/27678, 87308/27791, 87663/27904,
88018/28017, 88373/28130, 88728/28243, 89083/28356, 89438/28469, 89793/28582,
90148/28695, 90503/28808, 90858/28921, 91213/29034, 91568/29147, 91923/29260,
92278/29373, 92633/29486, 92988/29599, 93343/29712, 93698/29825, 94053/29938,
94408/30051, 94763/30164, 95118/30277, 95473/30390, 95828/30503, 96183/30616,
96538/30729, 96893/30842, 97248/30955, 97603/31068, 97958/31181, 98313/31294,
98668/31407, 99023/31520, 99378/31633, 99733/31746, 100088/31859, 100443/31972,
100798/32085, 101153/32198, 101508/32311, 101863/32424, 102218/32537, 102573/32650,
102928/32763, 103283/32876, 103638/32989, 103993/33102, 104348/33215, 208341/66317,
312689/99532, 833719/265381, 1146408/364913, 3126535/995207,
4272943/1360120, 5419351/1725033, 42208400/13435351, 47627751/15160384,
53047102/16885417, 58466453/18610450, 63885804/20335483, 69305155/22060516,
74724506/23785549, 80143857/25510582,
所有这些术语都是正确的，尽管如果增加最大值，则由于精度而开始遇到错误。我本人对只有13个汇合的术语给我留下了深刻的印象。 （如您所见，有一个小错误，有时它不会打印出第一个“ N/1”近似值，或者不正确地打印出来 - 我将其留给您解决！）
您可以尝试使用√2，其持续分数为[1; 2、2、2、2…]：
% ./a.out 1.41421356237309504880
1.414213562373095
1: 1/1
2: 3/2
2: 7/5
2: 17/12
2: 41/29
2: 99/70
2: 239/169
2: 577/408
2: 1393/985
2: 3363/2378
2: 8119/5741
2: 19601/13860
2: 47321/33461
3/2, 4/3, 7/5, 17/12, 24/17, 41/29, 99/70, 140/99, 239/169, 577/408, 816/577, 1393/985, 3363/2378, 4756/3363, 8119/5741, 19601/13860, 47321/33461,
或黄金比φ=（1+√5）/2的持续分数为[1; 1、1、1，…]：
% ./a.out 1.61803398874989484820
1.618033988749895
1: 1/1
1: 2/1
1: 3/2
1: 5/3
1: 8/5
1: 13/8
1: 21/13
1: 34/21
1: 55/34
1: 89/55
1: 144/89
1: 233/144
1: 377/233
2/1, 3/2, 5/3, 8/5, 13/8, 21/13, 34/21, 55/34, 89/55, 144/89, 233/144, 377/233,
（请参阅fibonacci数字？在这里，收敛是所有近似值。）
或有理性数字，例如4/3 = [1; 3]：
% ./a.out 1.33333333333333333333
1.333333333333333
1: 1/1
3: 4/3
3/2, 4/3,
或14/11 = [1; 3、1、2]：
% ./a.out 1.27272727272727272727
1.272727272727273
1: 1/1
3: 4/3
1: 5/4
2: 14/11
3/2, 4/3, 5/4, 9/7, 14/11,
享受！
其他提示
除了您无法相信其余部分的检查外，C程序还不错，从计算X*PQ也可以看出：
Iteration #1: 3: 3/1 - delta: 0.141592653589793116, rem: 0.141592653589793116
Iteration #2: 7: 22/7 - delta: -0.008851424871448188, rem: 0.062513305931051878
Iteration #3: 15: 333/106 - delta: 0.008821280518070296, rem: 0.996594406684156776
Iteration #4: 1: 355/113 - delta: -0.000030144353377892, rem: 0.003417231014946418
Iteration #5: 292: 103993/33102 - delta: 0.000019129331725765, rem: 0.634590879621879211
Iteration #6: 1: 104348/33215 - delta: -0.000011015021655680, rem: 0.575818424298580694
Iteration #7: 1: 208341/66317 - delta: 0.000008114310077190, rem: 0.736658567704091524
Iteration #8: 1: 312689/99532 - delta: -0.000002900711592702, rem: 0.357480987585133375
Iteration #9: 2: 833719/265381 - delta: 0.000002312886920208, rem: 0.797351564778957706
Iteration #10: 1: 1146408/364913 - delta: -0.000000587824615650, rem: 0.254151925163927682
Iteration #11: 3: 4272943/1360120 - delta: 0.000000549413016415, rem: 0.934654436927838420
Iteration #12: 1: 5419351/1725033 - delta: -0.000000038411599235, rem: 0.069914142051204637
Iteration #13: 14: 80143857/25510582 - delta: 0.000000011648808140, rem: 0.303257833981669641
Iteration #14: 3: 245850922/78256779 - delta: -0.000000003463355824, rem: 0.297524047014214316
Iteration #15: 3: 817696623/260280919 - delta: 0.000000001280568540, rem: 0.361072861287829440
Iteration #16: 2: 1881244168/598818617 - delta: -0.000000000931322575, rem: 0.769524124392304913
Iteration #17: 1: 2698940791/859099536 - delta: 0.000000000232830644, rem: 0.299504418772708979
Iteration #18: 3: 9978066541/3176117225 - delta: 0.000000000000000000, rem: 0.338848902789946401 ******* 'true' deviation below epsilon threshold
不隶属于 StackOverflow
