算法挑战:为浮点数生成连续分数 [英] Algorithm Challenge: Generate Continued Fractions for a float
问题描述
(编辑:为回应脾气暴躁的评论，不，这不是功课.我正在研究音高检测，获取一系列潜在的谐波峰，并尝试构建用于基本频率，因此，这实际上是一个非常实际的问题.)
(EDIT: In response to grumpy comments, No it isn't homework. I am working on pitch detection, taking an array of potential harmonic peaks, and attempting to construct candidates for fundamental frequency. So, it is actually a very practical question.)
考虑pi的最佳分数近似，按分母递增的顺序排列:3/1，22/7，355/113，...
Consider the best fractional approximations for (eg) pi, ordered by increasing denominator: 3/1, 22/7, 355/113, ...
挑战:创建一个 tidy C算法，该算法将为给定的浮点数生成第n个商近似a/b，同时还返回差异.
Challenge: create a tidy C algorithm that will generate the n'th quotient approximation a/b for a given float, returning also the discrepancy.
calcBestFrac(float frac，int n，int * a，int * b，float * err){...}
calcBestFrac(float frac, int n, int * a, int * b, float * err) {...}
我认为最好的技术是连续分数
拿走pi的小数部分，您得到3
现在，余数为0.14159 ... = 1/7.06251 ..
Take away the fractional part of pi, and you get 3
Now, the remainder is 0.14159... = 1/7.06251..
因此，下一个最佳有理数是3 + 1/7 = 22/7
从7.06251中减去7，您将得到0.06251 ..大约是1/15.99659 ..
So the next best rational is 3 + 1/7 = 22/7
Take away the 7 from 7.06251 and you get 0.06251.. Roughly 1/15.99659..
将其称为16，那么下一个最佳逼近度是
3 + 1/(7 + 1/16)= 355/113
Call it 16, then the next best approximation is
3 + 1/(7 + 1/16) = 355/113
但是，将其转换为纯净的C代码绝非易事.如果我整理整齐，我会发布.同时，有人可能会喜欢脑筋急转弯.
However, this is far from trivial to convert into clean C code. I will post if I get something tidy. In the meanwhile, someone may enjoy it as a brainteaser.
推荐答案
[因为您要求的是答案而不是评论.]
对于任何实数，其连续分数的收敛点p [k]/q [k]始终是最佳有理逼近，但并不是所有的 最佳有理逼近.要获得所有这些，您还必须采用半收敛/中位数-对于某些n≥1的整数，其形式为(p[k]+n*p[k+1])/(q[k]+n*q[k+1])的分数.取n = a [k + 2]得出p [k + 2]/q [k + 2]，要取的整数n是下限(a [k + 2]/2)或上限(a [ k + 2]/2)到a [k + 2].在Wikipedia上 也已提及.
For any real number, the convergents p[k]/q[k] of its continued fraction are always best rational approximations, but they aren't all the best rational approximations. To get all of them, you also have to take the semi-convergents/mediants — fractions of the form (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) for some integer n≥1. Taking n=a[k+2] gives p[k+2]/q[k+2], and the integers n to take are those from either floor(a[k+2]/2) or ceiling(a[k+2]/2), to a[k+2]. This is also mentioned on Wikipedia.
π的继续分数是[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 ).
The continued fraction for π is [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…] (sequence A001203 in OEIS), the sequence of convergents is 3/1, 22/7, 333/106, 355/113, 103993/33102… (A002485/A002486), and the sequence of best approximations is 3/1, 13/4, 16/5, 19/6, 22/7, 179/57… (A063674/A063673).
因此，该算法表示π的最佳近似值. = [3; 7，15，1，292，1，1，...]是
So the algorithm says that the best approximations of π = [3; 7, 15, 1, 292, 1, 1,…] are
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找到连续的分数(将a []中的项以及p []和q []中的收敛–排除全局变量)，函数all_best显示所有最佳有理逼近. >
Program
Here's a C program that given a positive real number, generates its continued fraction, its convergents, and the sequence of best rational approximations. The function find_cf finds the continued fraction (putting the terms in a[] and the convergents in p[] and q[] — excuse the global variables), and the function all_best prints all the best rational approximations.
#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;
}
示例
对于&pi ;，这是该程序的输出，大约需要0.003秒(即，这比循环遍历所有可能的分母要好！)，并进行换行以提高可读性:
Examples
For π, here's the output of this program, in about 0.003 seconds (i.e., it's truly better than looping through all possible denominators!), line-wrapped for readability:
% ./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,
所有这些术语都是正确的，尽管如果增加MAX会因为精度而开始出现错误.我对只有13个聚合词获得多少个词感到印象深刻. (如您所见，有一个小错误，有时它不会打印出第一个``n/1''近似值，或者打印不正确-我留给您修复！)
All these terms are correct, though if you increase MAX you start getting errors because of precision. I'm myself impressed with how many terms you get with only 13 convergents. (As you can see, there's a small bug where it sometimes doesn't print the very first "n/1" approximation, or prints it incorrectly — I leave it for you to fix!)
您可以尝试使用√2，其连续分数为[1; 2，2，2，2…]:
You can try with √2, whose continued fraction is [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，…]:
Or the golden ratio φ = (1+√5)/2 whose continued fraction is [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,
(请参阅斐波那契数?这里收敛是所有的近似值.)
(See the Fibonacci numbers? Here the convergents are all the approximants.)
或者使用有理数，例如4/3 = [1; 3]:
Or with rational numbers like 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]:
or 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,
享受！
这篇关于算法挑战:为浮点数生成连续分数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
