Efficient floating point scaling in C++

Candidate improvements

• Convert the radians x to rotations by dividing by 2*pi.

• Retain only the fraction so we have an angle (-1.0 ... 1.0). This simplifies the OP's modulo step to a simple "drop the whole number" step instead. Going forward with different angle units simply involves a co-efficient set change. No need to scale back to radians.

• For positive values, subtract 0.5 so we have (-0.5 ... 0.5) and then flip the sign. This centers the possible values about 0.0 and makes for better convergence of the approximating polynomial as compared to the math sine function. For negative values - see below.

• Call my_sin_small1() that uses this (-0.5 ... 0.5) rotations range rather than [-pi ... +pi] radians.

• In my_sin_small1(), fold constants together to drop the corr * step.

• Rather than use the truncated Taylor's series, use a more optimal set. IMO, this will provide better answers, especially near +/-pi.

Notes: No int to/from float code. With more analysis, possible to get a better set of coefficients that fix my_sin(+/-pi) closer to 0.0. This is just a quick set of code to demo less FP steps and good potential results.

C like code for OP to port to C++

FLOAT my_sin_small1(FLOAT x) {
  static const FLOAT A1 = -5.64744881E+01;
  static const FLOAT A2 = +7.81017968E+01;
  static const FLOAT A3 = -4.11145353E+01;
  static const FLOAT A4 = +6.27923581E+00;
  const FLOAT x2 = x * x;
  return x * (x2 * (x2 * (x2 * A1 + A2) + A3) + A4);
}

FLOAT my_sin1(FLOAT x) {
  static const FLOAT pi = 3.141592653589793238462;
  static const FLOAT pi2i = 1/(pi * 2);
  x *= pi2i;
  FLOAT xfraction = 0.5f - (x - truncf(x));
  return my_sin_small1(xfraction);
}

For negative values, use -my_sin1(-x) or like code to flip the sign - or add 0.5 in the above minus 0.5 step.

Test

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

int main(void) {
  for (int d = 0; d <= 360; d += 20) {
    FLOAT x = d / 180.0 * M_PI;
    FLOAT y = my_sin1(x);
    printf("%12.6f %11.8f  %11.8f\n", x, sin(x), y);
  }
}

Output

0.000000  0.00000000  -0.00022483
0.349066  0.34202013   0.34221691
0.698132  0.64278759   0.64255589
1.047198  0.86602542   0.86590189
1.396263  0.98480775   0.98496443
1.745329  0.98480775   0.98501128
2.094395  0.86602537   0.86603642
2.443461  0.64278762   0.64260530
2.792527  0.34202022   0.34183803
3.141593 -0.00000009   0.00000000
3.490659 -0.34202016  -0.34183764
3.839724 -0.64278757  -0.64260519
4.188790 -0.86602546  -0.86603653
4.537856 -0.98480776  -0.98501128
4.886922 -0.98480776  -0.98496443
5.235988 -0.86602545  -0.86590189
5.585053 -0.64278773  -0.64255613
5.934119 -0.34202036  -0.34221727
6.283185  0.00000017  -0.00022483

Alternate code below makes for better results near 0.0, yet might cost a tad more time. OP seems more inclined to speed.

FLOAT xfraction = 0.5f - (x - truncf(x));

// vs.

FLOAT xfraction = x - truncf(x);
if (x >= 0.5f) x -= 1.0f;

[Edit]

Below is a better set with about 10% reduced error.

-56.0833765f
77.92947047f
-41.0936875f
6.278635918f

Yet another approach:

Spend more time (code) to reduce the range to ±pi/4 (±45 degrees), then possible to use only 3 or 2 terms of a polynomial that is like the usually Taylors series.

float sin_quick_small(float x) {
  const float x2 = x * x;
#if 0
  // max error about 7e-7
  static const FLOAT A2 = +0.00811656036940792f;
  static const FLOAT A3 = -0.166597759850666f;
  static const FLOAT A4 = +0.999994132743861f;
  return x * (x2 * (x2 * A2 + A3) + A4);
#else
  // max error about 0.00016
  static const FLOAT A3 = -0.160343346851626f;
  static const FLOAT A4 = +0.999031566686144f;
  return x * (x2 * A3 + A4);
#endif
}

float cos_quick_small(float x) {
  return cosf(x); // TBD code.
}

float sin_quick(float x) {
  if (x < 0.0) {
    return -sin_quick(-x);
  }
  int quo;
  float x90 = remquof(fabsf(x), 3.141592653589793238462f / 2, &quo);
  switch (quo % 4) {
    case 0:
      return sin_quick_small(x90);
    case 1:
      return cos_quick_small(x90);
    case 2:
      return sin_quick_small(-x90);
    case 3:
      return -cos_quick_small(x90);
  }
  return 0.0;
}

int main() {
  float max_x = 0.0;
  float max_error = 0.0;
  for (int d = -45; d <= 45; d += 1) {
    FLOAT x = d / 180.0 * M_PI;
    FLOAT y = sin_quick(x);
    double err = fabs(y - sin(x));
    if (err > max_error) {
      max_x = x;
      max_error = err;
    }
    printf("%12.6f %11.8f  %11.8f  err:%11.8f\n", x, sin(x), y, err);
  }
  printf("x:%.6f err:%.6f\n", max_x, max_error);
  return 0;
}