马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?欢迎注册
×
本文同步发表于我的csdn博客,见http://blog.csdn.net/liangbch/archive/2008/01/24/2064108.aspx
实数包含有理数和无理数,任何有理数都可以表示为p/q(p,q是整数,q>0)的形式,如果指定一个分数的分母不超过某个值,对于一般的有理数或者无理数,是不可以用一个分数来准确地表示的。我们这里主要讨论,如何找出一对分数p1/q1和p2/q2,使得q1 和q2 小于给定的值n,而p1/q1和p2/q2尽可能接近一个给定的实数.。为了便于说明,我们用C++语言的格式给出问题的定义:
函数接口void search( int &p1, int &q1, int &p2,int &q2, int n, double f)
功能: 找出一对不可约p1/q1 和p2/q2, 使得这两个分数的分母不大于n, 且p1/q1 <= f <= p2/q2, 且这两个分数尽可能接近f。
在解决这个问题之前,我们先介绍一些准备知识。
定义1:最简分数(也称既约分数或不可约分数)。若p,q的最大公约数是1,我们称分数p/q 是最简分数。不特别说明,以下提到的分数的分子和分母均为非负整数。
定义2:真分数,若p,q是正整数,0 0
固(a+c)/(b+d)>a/b。同理可证 (a+c)/(b+d) < c/d
定义3:法雷数列
对任意给定的一个自然数n,将分母小于等于n的不可约的真分数按升序排列,并且在第一个分数之前加上0/1,在最后一个分数之后加上1/1,这个序列称为n级法雷数列,以Fn表示。如F5为:1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5.
法雷数列的构造:
法雷数列的构造可采用2分法,即如果 a/b, c/d (a/b法雷数列的性质:- 除了1级法雷数列外,所有的法雷数列都有奇数个元素,其中居于正中间的那个元素一定是1/2.
- 当n趋于正无穷时,n级法雷数列包含的元素的个数趋于3/(π*π) * n2 ≈ 0.30396355 * n2.
- n级法雷数列中,若相邻两个元素是a/b 和c/d (a/b
实数化分数方法
对于有理数0,我们可以用0/1表示;对于有理数x<0, 总可以表示为 –(p/q), 其中p>0,q>0;而对于所有大于等于1的正有理数,总可以表示为 n + p/q ( n, p, q为非负整数,q!=1, p
= a/b且f<=c/d,a/b称为实数f的下界,c/d称为实数f的上界,求这个下界和上界实际上是找出一个n级法雷数列中两个相邻的元素。下面是化一个小于1的正实数为分数的算法。
Step1: 置实数f 的下界为 a/b=0/1, 上界为c/d =1/1。
Step2: 计算出下界和上界之间的数 p/q = (a+c)/(b+d)
Step3: 若 q>n(分母大于指定值),计算中止。
若 p/q =f,置a/b为p/q , 置c/d 为p/q, 计算中止。
若 p/q > f, 置下界a/b为p/q
若 p/q< f, 置上界c/d为p/q
Step4, 重复step 2-3
当计算终止时,a/b为这个实数的下界,c/d为这个实数的上界。
如果要转化的实数f小于1/2, 用上述逐步求精的步骤,计算出上下界,迭代次数稍多。我们可用下面的步骤代替step1, 直接找出一个更精确的上下界。 若e= 1/x 的向下取整,则这个实数的下界和上界为 1/(e+1)和1/e
误差分析
根据法雷数列性质3我们知道,n级法雷数列中相邻的两个元素可以表示一个区间 [a/b, c/d],前一个元素q/b为区间的下界,后一个元素c/d为区间的上界,这个区间的宽度h =c/d- a/b,满足 1/n <= h <1/(n*n-1)。若运气好的话,一个实数正好落在一个宽度为1/n(n-1) 的区间,这个区间的下界或上界与这个实数的差不超过abs(1/(n*(n-1)))。若运气很差,一个实数恰好小于法雷序列的第2个元素或者最后一个元素。则这个元素的下界和上界与这个实数的差不超过1/n。
下面给出一段示例代码,这段代码可计算出sqrt(2),sqrt(3),sqrt(5) 以及无理数π和e的 分数表示。
- #include "math.h"
- #include "stdlib.h"
- #include "stdio.h"
-
- typedef struct _frac
- {
- unsigned long numerator;
- unsigned long denominator;
- }FRAC;
-
- double getValue( FRAC f)
- {
- double t=(double)f.numerator / (double)f.denominator;
- return t;
- }
-
- //f必须为正数
- void searchFrac( FRAC* pLow, FRAC* pHigh,int n, double f)
- {
- FRAC low;
- FRAC high;
- FRAC mid;
-
- int k=(int)f;
- f -= (double)k;
-
- low.numerator=0;
- low.denominator=1;
- high.numerator=1;
- high.denominator=1;
-
- mid.numerator= low.numerator + high.numerator;
- mid.denominator=low.denominator + high.denominator;
-
- while ( mid.denominator < n && fabs(f-getValue(mid))>1e-15 )
- {
- if ( getValue(mid) > f)
- {
- high=mid;
- }
- else
- {
- low=mid;
- }
- mid.numerator= low.numerator + high.numerator;
- mid.denominator=low.denominator + high.denominator;
- }
-
- if (k>0)
- {
- low.numerator += k * low.denominator;
- high.numerator += k * high.denominator;
- }
- *pHigh = high;
- *pLow =low;
- }
-
- int main(int argc, char* argv[])
- {
- FRAC low, high;
- double t;
- double f[]=
- {
- 2,
- 3,
- 5,
- 3.1415926535897932384626433832795,
- 2.7182818284590452353602874713527
- };
-
- for (int i=0;i<5;i++)
- {
- if (i<3)
- t=sqrt(f[i]);
- else
- t=f[i];
-
- searchFrac(&low,&high,1000000,t);
- printf("f=%.16f,low=%u/%u,high=%u/%u,\tlow=%.15f high=%.15f\n",t,
- low.numerator,low.denominator,high.numerator,high.denominator,
- getValue(low),getValue(high));
- }
- return 0;
- }
复制代码
版权所有 liangbch@263.net, 转载或者引用请注明出处。
参考文献:
[ 本帖最后由 mathe 于 2008-1-25 08:27 编辑 ]
|