c - sin(pi)没有返回精确的零？

This question already has answers here:

Closed 2 years ago.

Is floating point math broken?
（31个答案）

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

int main(int argc, char* argv[])
{
    double pi = acos(-1);

    for(int i = 0; i<65;i++)
    {
        double resultD = sin( pi * (2.0 *i / 128.0 ));
        float resultF  = sin( pi * (2.0f*i / 128.0f));

        printf("sine[%d]\n %.25f\n %.25f \n", i,resultD,resultF);
    }
}

在上面的代码中，我希望最后一个值是干净的零。然而，它类似于0.0000000000000001224606354。
另一方面，i=0处的值为clean0.000000000000000，i=32处的值为clean1.000000000000000。
我正在将浮点代码移植到定点代码，这一小的差异会升级。
有没有办法强迫i=64处的值精确为零？

最佳答案

sin(pi)不返回0，因为pi不是π。
π是一个无理数。即使具有扩展精度和最好的acos()/sin()例程，π的值也只能近似。从acos(-1)返回的值是机器pi，最接近π的double。

printf("machine pi %.20f...\n", pi);
// machine pi 3.14159265358979311600...
// math pi    3.14159265358979323846...
//                             ^^^^^...

当代码计算pi * (2.0*64/128.0)时，可以预期产品是精确的，没有舍入误差，pi的值是机器pi，而不是π。机器π的正弦确实是1.22…e-16
我正在将浮点代码移植到定点代码，这一小的差异会升级。
对于三角函数，如果原始参数是度数（或在OP的情况下是圆的1/128）或某个精确周期，与无理的2*π弧度不同：首先以原始单位进行范围缩减，然后可能使用三角恒等式。
最后转换成弧度并应用trig函数。

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

#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795
#endif

// Binary angle measurement
#define BAM 128

static double bam128tor(int bam) {
  return M_PI * 2 / BAM * bam;
}

double sin_BAM128(int bam) {
  // mod by the integral number of units in a circle
  bam %= BAM;
  if (bam < 0) bam += BAM;

  unsigned quadrant = 0;
  if (bam >= BAM / 2) {
    bam -= BAM / 2;
    quadrant += 2;
  }
  if (bam >= BAM / 4) {
    bam -= BAM / 4;
    quadrant += 1;
  }
  if (bam >= BAM / 8) {
    bam = -(BAM / 4 - bam);
    quadrant = (quadrant + 1) % 4;
  }
  printf("%2u %4d ", quadrant, bam);
  switch (quadrant % 4) {
    case 0:
      return sin(bam128tor(bam));
    case 1:
      return cos(bam128tor(bam));
    case 2:
      return sin(bam128tor(-bam)) * 1.0;
    default:
      return -cos(bam128tor(bam));
  }
}

int main(void) {
  double pi = acos(-1);
  for (int i = 0; i <= BAM; i += 4) {
    double resultD = sin(pi * (2.0 * i / 128.0));
    //float resultF = sin(pi * (2.0f * i / 128.0f));
    printf("sine[%3d]\t% .20f\t% .20f\n", i, resultD, sin_BAM128(i));
  }
}

输出

 0    0 sine[  0]    0.00000000000000000000  0.00000000000000000000
 0    4 sine[  4]    0.19509032201612824808  0.19509032201612824808
 0    8 sine[  8]    0.38268343236508978178  0.38268343236508978178
 0   12 sine[ 12]    0.55557023301960217765  0.55557023301960217765
 1  -16 sine[ 16]    0.70710678118654746172  0.70710678118654757274
 1  -12 sine[ 20]    0.83146961230254523567  0.83146961230254523567
 1   -8 sine[ 24]    0.92387953251128673848  0.92387953251128673848
 1   -4 sine[ 28]    0.98078528040323043058  0.98078528040323043058
 1    0 sine[ 32]    1.00000000000000000000  1.00000000000000000000
 1    4 sine[ 36]    0.98078528040323043058  0.98078528040323043058
 1    8 sine[ 40]    0.92387953251128673848  0.92387953251128673848
 1   12 sine[ 44]    0.83146961230254545772  0.83146961230254523567
 2  -16 sine[ 48]    0.70710678118654757274  0.70710678118654746172
 2  -12 sine[ 52]    0.55557023301960206663  0.55557023301960217765
 2   -8 sine[ 56]    0.38268343236508989280  0.38268343236508978178
 2   -4 sine[ 60]    0.19509032201612860891  0.19509032201612824808
 2    0 sine[ 64]    0.00000000000000012246  0.00000000000000000000 exact
 2    4 sine[ 68]   -0.19509032201612835911 -0.19509032201612824808
 2    8 sine[ 72]   -0.38268343236508967076 -0.38268343236508978178
 2   12 sine[ 76]   -0.55557023301960195560 -0.55557023301960217765
 3  -16 sine[ 80]   -0.70710678118654746172 -0.70710678118654757274
 3  -12 sine[ 84]   -0.83146961230254523567 -0.83146961230254523567
 3   -8 sine[ 88]   -0.92387953251128651644 -0.92387953251128673848
 3   -4 sine[ 92]   -0.98078528040323031956 -0.98078528040323043058
 3    0 sine[ 96]   -1.00000000000000000000 -1.00000000000000000000
 3    4 sine[100]   -0.98078528040323043058 -0.98078528040323043058
 3    8 sine[104]   -0.92387953251128662746 -0.92387953251128673848
 3   12 sine[108]   -0.83146961230254545772 -0.83146961230254523567
 0  -16 sine[112]   -0.70710678118654768376 -0.70710678118654746172
 0  -12 sine[116]   -0.55557023301960217765 -0.55557023301960217765
 0   -8 sine[120]   -0.38268343236509039240 -0.38268343236508978178
 0   -4 sine[124]   -0.19509032201612871993 -0.19509032201612824808
 0    0 sine[128]   -0.00000000000000024493  0.00000000000000000000 exact