Mise à l'échelle efficace de la virgule flottante en C++

Question

Je travaille sur mon implémentation rapide (et précise) du péché en C++, et j'ai un problème concernant la mise à l'échelle efficace de l'angle dans la fonction +- pi/2 gamme.

Ma fonction sin pour +-pi/2 en utilisant les séries de Taylor est la suivante (Note : FLOAT est une macro étendue à float ou double juste pour le repère)

/**
 * Sin for 'small' angles, accurate on [-pi/2, pi/2], fairly accurate on [-pi, pi]
 */
// To switch between float and double
#define FLOAT float

FLOAT
my_sin_small(FLOAT x)
{
    constexpr FLOAT C1 = 1. / (7. * 6. * 5. * 4. * 3. * 2.);
    constexpr FLOAT C2 = -1. / (5. * 4. * 3. * 2.);
    constexpr FLOAT C3 = 1. / (3. * 2.);
    constexpr FLOAT C4 = -1.;
    // Correction for sin(pi/2) = 1, due to the ignored taylor terms
    constexpr FLOAT corr = -1. / 0.9998431013994987;

    const FLOAT x2 = x * x;

    return corr * x * (x2 * (x2 * (x2 * C1 + C2) + C3) + C4);
}

Jusqu'ici tout va bien... Le problème survient lorsque j'essaie de mettre à l'échelle un angle arbitraire dans la plage +-pi/2. Ma solution actuelle est la suivante :

FLOAT
my_sin(FLOAT x)
{
    constexpr FLOAT pi = 3.141592653589793238462;
    constexpr FLOAT rpi = 1 / pi;

    // convert to +-pi/2 range
    int n = std::nearbyint(x * rpi);

    FLOAT xbar = (n * pi - x) * (2 * (n & 1) - 1);
    // (2 * (n % 2) - 1) is a sign correction (see below)
    return my_sin_small(xbar);
};

J'ai fait un repère et je perds beaucoup pour la mise à l'échelle +-pi/2.
Tricking avec int(angle/pi + 0.5) est un nope puisque c'est limité à la précision de l'int, nécessite également +- branchements, et j'essaie d'éviter les branches ...
Que dois-je essayer pour améliorer les performances de cette mise à l'échelle ? Je suis à court d'idées.

Résultats des tests de référence pour le flotteur. (Dans le benchmark, l'angle pourrait être hors de la plage de validité de my_sin_small, mais pour le bench, je ne m'en soucie pas...) :

Résultats du benchmarking pour le double.

Correction du signe pour xbar dans my_sin() :

Précision de l'algorithme par rapport à la fonction python sin() :

Answer 1

Améliorations apportées par les candidats

• Convertir les radians x a rotations en divisant par 2*pi.

• Ne retenez que la fraction pour obtenir un angle (-1,0 ... 1,0). Cela simplifie le problème de l'OP modulo à une simple étape de "suppression du nombre entier". Aller de l'avant avec différentes unités d'angle implique simplement un changement de jeu de coefficients. Il n'est pas nécessaire de revenir aux radians.

• Pour les valeurs positives, soustrayez 0,5 de manière à obtenir (-0,5 ... 0,5), puis inversez le signe. Cela permet de centrer les valeurs possibles autour de 0,0 et d'améliorer la convergence du polynôme d'approximation par rapport à l'algorithme sine fonction. Pour les valeurs négatives - voir ci-dessous.

• Appelez my_sin_small1() qui utilise cette plage de rotations (-0,5 ... 0,5) plutôt que [-pi ... +pi] radians.

• Sur my_sin_small1() , pliez les constantes ensemble pour faire tomber le corr * étape.

• Plutôt que d'utiliser la série de Taylor tronquée, utilisez un ensemble plus optimal. IMO, cela fournira de meilleures réponses, en particulier près de +/-pi.

Notes : Non int vers/depuis float code. Avec plus d'analyse, il est possible d'obtenir un meilleur ensemble de coefficients qui fixent my_sin(+/-pi) plus proche de 0,0. Il s'agit juste d'un ensemble de code rapide pour démontrer les étapes moins FP et les bons résultats potentiels.

---

Code C pour que l'OP puisse le porter en 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);
}

Pour les valeurs négatives, utilisez -my_sin1(-x) ou un code similaire pour inverser le signe - ou ajouter 0,5 dans l'étape moins 0,5 ci-dessus.

---

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

Sortie

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

---

Le code alternatif ci-dessous permet d'obtenir de meilleurs résultats près de 0.0, mais il peut coûter un peu plus de temps. L'OP semble plus enclin à la rapidité.

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

// vs.

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

---

[Edit]

Vous trouverez ci-dessous un meilleur ensemble avec une erreur réduite d'environ 10%.

-56.0833765f
77.92947047f
-41.0936875f
6.278635918f

Answer 2

Une autre approche encore :

Passez plus de temps (code) pour réduire la gamme à ±pi/4 (±45 degrés), alors possible d'utiliser seulement 3 ou 2 termes d'un polynôme qui est comme la série de Taylors habituellement.

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

Answer 3

   FLOAT
   my_sin(FLOAT x)
   {
       constexpr FLOAT pi = 3.141592653589793238462;
       constexpr FLOAT rpi = 1 / pi;

       // convert to +-pi/2 range
       int n = std::nearbyint(x * rpi);

       FLOAT xbar = ((float)n * pi - x) * (2.0 * (n % 2.0) - 1.0);
       // (2 * (n % 2) - 1) is a sign correction (see below)
       return my_sin_small(xbar);
   };

Essayez 2.0 au lieu de 2