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

Question

这是我需要编程的数学计算 - 将序列 a_n 定义为

• a₁ = 1
• a_(n+1) = 1 / (2 × [a_n] - a_n + 1) 其中 [x]是floor 函数

如果a₂₀₁₄ = 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;
}

问题；

• 我无法编译（使用网页编译器，例如键盘等）

• 不知道怎么出结果

• 分数结果错误

Answer 1

这是一个危险的序列：

floor(具有潜在表示错误的值) = 大错误！！

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) 替换为floor(x + 2*expectedError)

2) 关于累积错误：将double 表示替换为(v1:int, v2:int) 分数表示。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
}

\谢谢{SinanÜnür}

Answer 2

对于总是喜欢提醒人们注意浮点错误的人，在这种情况下，我在方向盘上睡着了。

如果你运行以下程序，

/* 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

但是，这是由于 @JJoao points 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 语言提供了各种工具来操作浮点数的四舍五入方式。在这种特殊情况下，鉴于 floor 函数的性质，将舍入模式设置为 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