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

Question

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.

Answer 1

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

Answer 2

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

This part is somewhat easy. The approximate error in limiting sine(x) code to N polynomial terms is about the value of the N+1^th term of the Taylor series.

For a limited range of x [0 ... π/4] about 3 terms needed for 15 bits. See sine series definitions.

4th term is (π/4)⁷/7! or 3.7e-5 or about 1 part in 2^14.7. Trig identities can handle the rest of the number range.

---

In trying to create a uint16_t sin_f(uint16_t) an early problem is how to scale the input and output.

The result of math sine is [-1.0 ... +1.0]. By only calculating the first 90 degrees on sine, the range is [0.0 ... +1.0]. Scaling this by a power-of-2 could use [0 ... 65536], but the ends are inclusive so that result will not fit in a uint16_t. Perhaps use [0 ... 32768]?

OP implied 1.0 is one revolution, 360 degrees or 2*π radians. So input is a 16-bit fraction [0 ... 65535].

---

Below is a modest attempt that uses a 4 term polynomial. Terms were found via some excel curve fitting techniques and are not necessarily optimal. Its has 2 known problems: result in the range [0 ... 32769] (could tweak the scale a bit to fix) and worst case off-by-4 result. (off by 2 was my goal.) It does offer some idea to OP of what is involved. As other say, this is not trivial and a dynamic solution to variant fixed widths looks to be very hard. A constant 16 bit fixed point was hard enough.

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>

uint16_t mul16(uint16_t a, uint16_t b, int shift) {
  uint32_t y32 = a;
  y32 *= b;
  y32 >>= shift - 1;
  if (y32 & 1) {
    y32 >>= 1;
    y32 += 1;
  } else {
    y32 >>= 1;
  }

  if (y32 > 0xFFFF) {
    printf("@@@@@@ %u %u %llu %d\n", 1u*a, 1u*b, 1llu*y32, shift);
    exit(0);
  }
  uint16_t y16 = (uint16_t) y32;
  return y16;
}

// 0 to 0x4000  (map to 0 to 90 degrees or 0 to pi/2 R)
// pseudo code:  = x*(1 - b*x^2 + c*x^4 - d*x^6)
uint16_t sine_fixed(uint16_t x) {
  const uint16_t t3 = 53902; // 431214.77/8;
  const uint16_t t5 = 64636; // 129272.18/2;
  const uint16_t t7 = 58833; // 58833.22

  uint16_t xx = mul16(x, x, 16);
  uint16_t term = 51472; // 2*pi*65536 / 8
  uint16_t sum = term;

  term = mul16(mul16(term, xx, 16 - 3), t3, 16 - 0);
  sum = (uint16_t) (sum - term);

  term = mul16(mul16(term, xx, 16 - 1), t5, 16 - 0);
  sum = (uint16_t) (sum + term);

  term = mul16(mul16(term, xx, 16 - 0), t7, 16 - 0);
  sum = (uint16_t) (sum - term);

  uint16_t y = mul16(x, sum, 16 - 3 + 1);
  return y;
}

Test code

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

void sin_fixed_test(uint16_t x) {
  double pi = acos(-1);
  uint16_t y = sine_fixed(x);
  double X = x * 2*pi / 65536.0;
  double Y = sin(X);
  long ye =  lrint(Y * 65536.0/2);
  //printf("sine(%5u) --> %5u, expect %5ld\n", 1u * x, 1u * y, ye);
  long diff = labs(ye - y);
  static long diff_max = 0;
  if (diff > diff_max) {
    diff_max = diff;
    printf("sine(%5u) --> %5u, expect %5ld !!!\n", 1u * x, 1u * y, ye);
  }
}

void sin_fixed_tests() {
  sin_fixed_test(15887);
  for (uint16_t x = 0; x <= 65536u / 4u; x += 1) {
    sin_fixed_test(x);
  }
}

int main() {
  sin_fixed_tests();
  return 0;
}