为什么涉及底线函数的这个公式没有给出我期望的结果？

Question

这是数学运算，我需要程序 - 定义序列的_Ñ作为

• 一个 = 1
• 一个_{（N + 1）} = 1 /（2 × [一_ñ] - 一个_ñ + 1）其中，[X]是floor
功能

如果 = P/Q，发现p + q

而且，这是我尝试的代码：

#include <stdio.h> 
#include <math.h> 

main() { 
    int n=1; 
    double an=1.0; 
    double an_floor=0.0; 

    while (n<=2014) { 
     an_floor = floor(an); 
     an = 1/(2*an_floor-an+1); 
     n = n + 1; 
    } 

    printf("%lf", an); 

    return 0; 
} 

问题;

• （使用Web编者如键盘等），我不能编译

• 不知道如何使结果

• 分数结果是错误的

Answer 1

对于那些总是喜欢提醒人们浮点错误的人，在这种情况下，我被抓住了。

如果您运行下面的程序，

/* a1 = 1, a(n+1) = 1/(2*[an]-an+1) ([x] is floor function) 
* a2014 = p/q 
* find p+q 
*/ 

#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 

#define LAST_ELEMENT 2014 

static double 
next_element(double prev) { 
    return 1.0/(2.0 * floor(prev) - prev + 1.0); 
} 


int main(void) { 
    int i = 0; 
    double last = 1.0; 

    for (i = 0; i < LAST_ELEMENT; i += 1) { 
     last = next_element(last); 
    } 

    printf("%g\n", last); 
    return EXIT_SUCCESS; 
} 

你得到这样的输出：

C:\...\Temp> cl /fp:precise /O2 rrr.c 
C:\...\Temp> rrr.exe 
2.25

但是，这是由于浮点误差@JJoaopoints out。他概述了可以处理这个特定问题的具体方法。

另一种方法是利用任意精度的数学库来帮助你。首先，让我们用一个快速的Perl脚本验证问题：

#!/usr/bin/env perl 

use strict; 
use warnings; 

use POSIX qw(floor); 

sub next_element { '1'/(('2' * floor($_[0])) - $_[0] + '1.0') } 

sub main { 
    my ($x, $n) = @_; 
    $x = next_element($x) for 1 .. $n; 
    printf("%g\n", $x); 
} 

main(1, 2014); 

输出：

C:\...\Temp> perl t.pl 
2.25

同C程序。

现在，使用Perl的bignum：

C:\...\Temp> perl -Mbignum t.pl 
5.83333

这更高的精度是有性能开销。如果没有bignum，脚本运行时间为0.125秒，其中约0.094用于计算之外的其他事情。用bignum，大概需要两秒钟。

现在，现代C提供了各种设施来操纵浮动数字的四舍五入。在这种特殊情况下，考虑到地板函数的性质，舍入模式设置为FE_DOWNWARD将解决这个问题：

#include <fenv.h> 
#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 

#define LAST_ELEMENT 2014 

static double 
next_element(double prev) { 
    return 1.0/(2.0 * floor(prev) - prev + 1.0); 
} 


int main(void) { 
    int i = 0; 
    double last = 1.0; 

    const int original_rounding = fegetround(); 
    fesetround(FE_DOWNWARD); 

    for (i = 0; i < LAST_ELEMENT; i += 1) { 
     last = next_element(last); 
    } 

    fesetround(original_rounding); 
    printf("%g\n", last); 

    return EXIT_SUCCESS; 
} 

输出：

C:\...\Temp> cl /O2 /fp:precise rrr.c    
/out:rrr.exe               

C:\...\Temp> rrr         
5.83333

Answer 2

这是一个危险的序列：

地板（具有潜在代表性误差的值）=大错误！

floor(3-epsilon)= 2 
floor(3+epsilon)= 3 

大数累积实数计算=大错误!!

您需要研究数值精度错误！使用

的OP代码，我们得到

1 
0.5 
2 
0.333333 
1.5 
0.666667 
3 
0.5   <===== Big precision ERROR it should be 0.25 
2 
1 
2.2518e+15A <===== infinity 
4.44089e-16 <===== zero 
1 
0.5 
....(wrong) cyclic behavior 

最终（错误！）结果= 2.2518e+15A <===== infinity

更新：但是，我们该怎么办？

1）floor问题：关于累积误差由floor(x + 2*expectedError)

2）代替floor(x)：由(v1:int, v2:int) 分数表示替换double表示。 An -> v1/v2

newAn = newv1/newv2 = 
     = 1/(2 * flo(An) + 1 - An) = 
     = 1/(2 * flo(v1/v2) + 1 - v1/v2) = 
     = ... = 
     = v2/(2 * v2 * flo(v1/v2) + v2 - v1) 

这样

newv1 = v2 
newv2 = 2 * v2 * flo(v1/v2) + v2 - v1 

在C：

main() { 
    int v1=1, v2=1, n=1, oldv1; 
    double an_floor; 

    while (n<=2014) { 
    an_floor = floor((0.0+v1)/v2 + 0.0000000001); 

    oldv1=v1; 
    v1=v2; 
    v2=an_floor * v2 * 2 + v2 - oldv1; 
    n++; 
    // printf("%g\n", (0.0+v1)/v2); 
    } 

    printf("%d/%d=%g\n", v1, v2, (0.0+v1)/v2); //  35/6=5.83333 
} 

\ {感谢思南Ünür}