Calculating sine to full precision of given fixed point type, when to stop?

I'm writing a sine implementation for my own fixed-point library without using pre-computed tables or library functions, so just basic mathematical operations (add, subtract, multiply, divide). I want to calculate this sine value to the full precision of whatever fixed point type the result is in, e.g. one with 16 fractional bits.

As far as I know the Taylor series is the only way of doing this, and it gets more precise the more terms I add, but how do I determine when to stop adding terms? Is it enough to just check if the next term would be smaller than the target precision? Is it even practical to use the Taylor series in this way or do I need to use something else?

I'm using C and want to make the number of fractional bits of my fixed-point type (or types) configurable, which is why I need to be able to generalize my stopping condition in this way.

1 answer

  • answered 2018-02-02 18:44 chux

    Following is not a fixed point solution, but a floating point one. So it may provide insight for a fixed point one.

    It does use a library function remainder() for range reduction, but that could also be replaced once detailed coding goals are expressed.

    It uses recursion to add the smaller terms together first. With higher precision floating point, this may make for a deep stack excursion.

    The recursion termination here is a test with 1.0, which makes sense for floating point. I'd expect a compare against an epsilon for fixed point.

    static double sin_r_helper(double term, double xx, unsigned i) {
      if (1.0 + term == 1.0) return term;
      return term - sin_r_helper(term * xx / ((i + 1) * (i + 2)), xx, i + 2);
    }
    
    double sin_r(double x /* radians */) {
      x = remainder(x, 2 * M_PI);
      return x * sin_r_helper(1.0, x * x, 1);
    }