sin(pi) not returning exact zero?

Question

#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);
    }
}

In the above code, I would expect to get last value to be clean zero. However, it is something like 0.0000000000000001224606354. On the other hand, the values at i=0 is clean 0.000000000000000, and the value at i=32 is clean 1.000000000000000.

I am porting floating-point to fixed-point code and this small difference escalates.
Is there a way to force the value at i=64 to be exact zero?

And a precision of 1e-15 isn't good enough? It **can't** be absolutely precise in every case, see linked duplicate. — , Oct 10 '17 at 07:09
If you're performing large numbers of chained calculations, then I suspect that somehow constraining your `sin(pi)` result to zero will increase, rather than decrease, the total error. What happens if you take the sine of a number that is just a tiny bit on either side of `pi`? The different between this and the sine of `pi` will be larger than it ought to be. Get an abitrary-precision math library, and you can have a precision of a hundred digits, or a thousand, or whatever, at your discretion. — Kevin Boone, Oct 10 '17 at 07:23
@Danijel but the answer remains the same: this is how floating-point arithmetics work. I don't understand this, you're calculating `sin` with some custom-defined decimal format? — , Oct 10 '17 at 07:24
In my fixed-point implementation I have a table where `sin[64]=0`, and in floating-point version, like above, I get something that is not zero. This small difference further propagates through the code, and eventualy makes the final results of fixed-point implementation much different than floating-point. — Danijel, Oct 10 '17 at 07:27
If you plan to do chained floating point math, you need to analyze the effect of error propagation and decide on a suitable strategy for your case. Some times, a calculation can be made more robust just by changing the order of computing partial results (robust = decreasing the effect of computation errors instead of escalating it). Some times, the result is known to be a limited precision number and rounding sub-results to limited precision can lead to better results. Some times, you can only get better results by inreasing the precision and `float` is just the wrong type for the job. — grek40, Oct 10 '17 at 07:53
Well, again, this is how floating-point works. I actually don't see how trigonometric functions in fixed-point would be useful. Fixed-point is typically used when your data actually has a fixed decimal precision and your operations stay within that (like, when values representing money). Calculating trigonometrics will probably give better overall results with floating-point. — , Oct 10 '17 at 07:54
Right, but the processor does not have a Floating Point unit. — Danijel, Oct 10 '17 at 07:57
I smell an [XY Problem](https://meta.stackexchange.com/questions/66377/what-is-the-xy-problem)... how about a question that focuses on your actual issues with a non-fpu platform? — grek40, Oct 10 '17 at 08:00
See [What Every Computer Scientist should know about Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) — Chris Dodd, Oct 10 '17 at 14:16
this isn't a correct duplicate here, because the real problem is that **there's no floating-point value exact to pi**. Some candidates [`Sin` and `Cos` Functions returning incorrect results](https://stackoverflow.com/q/13846079/995714), [How can I work around the fact that in C++, `sin(M_PI)` is not 0?](https://stackoverflow.com/q/2670999/995714), [sin, cos, tan and rounding error](https://stackoverflow.com/q/1527588/995714) — phuclv, Oct 10 '17 at 14:54
@ChrisDodd Not a proper dupe. This issue here is specific: trigonometric argument range reduction and not general issue about Floating-point. This issue applies to OP's fixed point concern's too. As OP's original angle measure are exact with an exact period, high quality answers are possible. — chux - Reinstate Monica, Oct 10 '17 at 16:19

chux - Reinstate Monica · Answer 1 · 2017-10-10T15:37:39.347

sin(pi) does not return 0 because pi is not π.

π is an irrational number. Even with extended precision and the best acos()/sin() routines, the value of π can only be approximated. The value returned from acos(-1) is machine-pi, the closest double to π.

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

When code calculates pi * (2.0*64/128.0), the product can be expected to be exact with no rounding error, yet the value of pi here is machine-pi and not π. The sine of machine-pi is indeed about 1.22...e-16

I am porting floating-point to fixed-point code and this small difference escalates.

For trigonometric functions, if the original argument is degrees (or in OP's case 1/128th of a circle) or some exact period, unlike the irrational 2*π radians: perform the range reductions first in original units, then perhaps then use trig identities. Lastly convert to radians and apply the trig function.

#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));
  }
}

Output

 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

sin(pi) not returning exact zero?

1 Answers1

Related