TA貢獻(xiàn)1865條經(jīng)驗(yàn) 獲得超7個(gè)贊
我發(fā)現(xiàn)David Eppstein 發(fā)現(xiàn)給定實(shí)數(shù) C代碼的理性近似正是你所要求的。它基于連續(xù)分?jǐn)?shù)理論，非?？焖偾蚁喈?dāng)緊湊。
我已經(jīng)為特定的分子和分母限制使用了這個(gè)版本。
/*** find rational approximation to given real number** David Eppstein / UC Irvine / 8 Aug 1993**** With corrections from Arno Formella, May 2008**** usage: a.out r d** r is real number to approx** d is the maximum denominator allowed**** based on the theory of continued fractions** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))** then best approximation is found by truncating this series** (with some adjustments in the last term).**** Note the fraction can be recovered as the first column of the matrix** ( a1 1 ) ( a2 1 ) ( a3 1 ) ...** ( 1 0 ) ( 1 0 ) ( 1 0 )** Instead of keeping the sequence of continued fraction terms,** we just keep the last partial product of these matrices.*/#include <stdio.h>main(ac, av)int ac;char ** av;{
double atof();
int atoi();
void exit();
long m[2][2];
double x, startx;
long maxden;
long ai;
/* read command line arguments */
if (ac != 3) {
fprintf(stderr, "usage: %s r d\n",av[0]); // AF: argument missing
exit(1);
}
startx = x = atof(av[1]);
maxden = atoi(av[2]);
/* initialize matrix */
m[0][0] = m[1][1] = 1;
m[0][1] = m[1][0] = 0;
/* loop finding terms until denom gets too big */
while (m[1][0] * ( ai = (long)x ) + m[1][1] <= maxden) {
long t;
t = m[0][0] * ai + m[0][1];
m[0][1] = m[0][0];
m[0][0] = t;
t = m[1][0] * ai + m[1][1];
m[1][1] = m[1][0];
m[1][0] = t;
if(x==(double)ai) break; // AF: division by zero
x = 1/(x - (double) ai);
if(x>(double)0x7FFFFFFF) break; // AF: representation failure
}
/* now remaining x is between 0 and 1/ai */
/* approx as either 0 or 1/m where m is max that will fit in maxden */
/* first try zero */
printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
startx - ((double) m[0][0] / (double) m[1][0]));
/* now try other possibility */
ai = (maxden - m[1][1]) / m[1][0];
m[0][0] = m[0][0] * ai + m[0][1];
m[1][0] = m[1][0] * ai + m[1][1];
printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
startx - ((double) m[0][0] / (double) m[1][0]));}
TA貢獻(xiàn)1799條經(jīng)驗(yàn) 獲得超8個(gè)贊
從Python 2.6開始就有fractions模塊。
（引用文檔。）
>>> from fractions import Fraction>>> Fraction('3.1415926535897932').limit_denominator(1000)Fraction(355, 113)>>> from math import pi, cos>>> Fraction.from_float(cos(pi/3))Fraction(4503599627370497, 9007199254740992)>>> Fraction.from_float(cos(pi/3)).limit_denominator()Fraction(1, 2)
TA貢獻(xiàn)1752條經(jīng)驗(yàn) 獲得超4個(gè)贊
如果輸出是為了給人類讀者一個(gè)結(jié)果順序的快速印象，那么返回“113/211”之類的東西是沒有意義的，所以輸出應(yīng)該限制為使用一位數(shù)字（也許是1 / 10和9/10）。如果是這樣，您可以觀察到只有27 種不同的分?jǐn)?shù)。
由于用于生成輸出的基礎(chǔ)數(shù)學(xué)將永遠(yuǎn)不會(huì)改變，因此解決方案可能是簡單地對(duì)二進(jìn)制搜索樹進(jìn)行硬編碼，以便該函數(shù)最多執(zhí)行l(wèi)og（27）?= 4 3/4比較。這是經(jīng)過測(cè)試的C版代碼
char *userTextForDouble(double d, char *rval){
if (d == 0.0)
return "0";
// TODO: negative numbers:if (d < 0.0)...
if (d >= 1.0)
sprintf(rval, "%.0f ", floor(d));
d = d-floor(d); // now only the fractional part is left
if (d == 0.0)
return rval;
if( d < 0.47 )
{
if( d < 0.25 )
{
if( d < 0.16 )
{
if( d < 0.12 ) // Note: fixed from .13
{
if( d < 0.11 )
strcat(rval, "1/10"); // .1
else
strcat(rval, "1/9"); // .1111....
}
else // d >= .12
{
if( d < 0.14 )
strcat(rval, "1/8"); // .125
else
strcat(rval, "1/7"); // .1428...
}
}
else // d >= .16
{
if( d < 0.19 )
{
strcat(rval, "1/6"); // .1666...
}
else // d > .19
{
if( d < 0.22 )
strcat(rval, "1/5"); // .2
else
strcat(rval, "2/9"); // .2222...
}
}
}
else // d >= .25
{
if( d < 0.37 ) // Note: fixed from .38
{
if( d < 0.28 ) // Note: fixed from .29
{
strcat(rval, "1/4"); // .25
}
else // d >=.28
{
if( d < 0.31 )
strcat(rval, "2/7"); // .2857...
else
strcat(rval, "1/3"); // .3333...
}
}
else // d >= .37
{
if( d < 0.42 ) // Note: fixed from .43
{
if( d < 0.40 )
strcat(rval, "3/8"); // .375
else
strcat(rval, "2/5"); // .4
}
else // d >= .42
{
if( d < 0.44 )
strcat(rval, "3/7"); // .4285...
else
strcat(rval, "4/9"); // .4444...
}
}
}
}
else
{
if( d < 0.71 )
{
if( d < 0.60 )
{
if( d < 0.55 ) // Note: fixed from .56
{
strcat(rval, "1/2"); // .5
}
else // d >= .55
{
if( d < 0.57 )
strcat(rval, "5/9"); // .5555...
else
strcat(rval, "4/7"); // .5714
}
}
else // d >= .6
{
if( d < 0.62 ) // Note: Fixed from .63
{
strcat(rval, "3/5"); // .6
}
else // d >= .62
{
if( d < 0.66 )
strcat(rval, "5/8"); // .625
else
strcat(rval, "2/3"); // .6666...
}
}
}
else
{
if( d < 0.80 )
{
if( d < 0.74 )
{
strcat(rval, "5/7"); // .7142...
}
else // d >= .74
{
if(d < 0.77 ) // Note: fixed from .78
strcat(rval, "3/4"); // .75
else
strcat(rval, "7/9"); // .7777...
}
}
else // d >= .8
{
if( d < 0.85 ) // Note: fixed from .86
{
if( d < 0.83 )
strcat(rval, "4/5"); // .8
else
strcat(rval, "5/6"); // .8333...
}
else // d >= .85
{
if( d < 0.87 ) // Note: fixed from .88
{
strcat(rval, "6/7"); // .8571
}
else // d >= .87
{
if( d < 0.88 ) // Note: fixed from .89
{
strcat(rval, "7/8"); // .875
}
else // d >= .88
{
if( d < 0.90 )
strcat(rval, "8/9"); // .8888...
else
strcat(rval, "9/10"); // .9
}
}
}
}
}
}
return rval;}舉報(bào)