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

Posted 2023-02-16

【中文标题】为啥这个涉及 floor 函数的公式没有给出我期望的结果？

【英文标题】：Why is this formula involving the floor function not giving the results I expect?为什么这个涉及 floor 函数的公式没有给出我期望的结果？

【发布时间】：2015-05-07 09:57:59

【问题描述】：

这是我需要编程的数学计算 - 将序列 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;

问题；

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

不知道怎么出结果

分数结果错误

【问题讨论】：

递归在这里可能是一个更好的选择。

仔细分析，这其实是一个极好的问题！

【参考方案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

【讨论】：

请查看我的答案以获得替代修复。

【参考方案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

【讨论】：

会是2.2518e+15 = infinty吗？

非常有趣！！！ 1) 在 unix 中我无法重现 C 版本的行为。 （我得到2.5 （我不知道/fp:precise 对应的选项是什么）

Hmmmm ...你可以尝试添加#pragma STDC FENV_ACCESS ON吗？

将FE_DOWNWARD 替换为FE_UPWARD 即可！奇怪的是，我们已经在结果精度中丢失了 5 个十进制数字...... 5.8333333333428480 5.83333

思南，看起来像“超额浮点精度”、“舍入”和类似的地方是一个困难的领域，（错误、bw 编译器可能存在的差异等）。在我们的问题中：（1）fesetround(FE_UPWARD) 修复了地板问题（优秀） - 并且可能略微增加了全局最终错误（非常小的问题）。 (2) perl -Mbignum 有效，并且是一种非常清晰和通用的技术（感谢您告诉我有关 big* perl 模块的信息）。 (3) perl -Mbigrat 有效，是一个无表示错误解决方案。