C: Comment envelopper un float dans l’intervalle [-pi, pi)

Je recherche un code C agréable qui accomplira efficacement:

while (deltaPhase >= M_PI) deltaPhase -= M_TWOPI; while (deltaPhase < -M_PI) deltaPhase += M_TWOPI; 

Quelles sont mes options?

Edit 19 avril 2013:

La fonction Modulo a été mise à jour pour gérer les cas limites comme indiqué par aka.nice et arr_sea:

 static const double _PI= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348; static const double _TWO_PI= 6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696; // Floating-point modulo // The result (the remainder) has same sign as the divisor. // Similar to matlab's mod(); Not similar to fmod() - Mod(-3,4)= 1 fmod(-3,4)= -3 template T Mod(T x, T y) { static_assert(!std::numeric_limits::is_exact , "Mod: floating-point type expected"); if (0. == y) return x; double m= x - y * floor(x/y); // handle boundary cases resulted from floating-point cut off: if (y > 0) // modulo range: [0..y) { if (m>=y) // Mod(-1e-16 , 360. ): m= 360. return 0; if (m<0 ) { if (y+m == y) return 0 ; // just in case... else return y+m; // Mod(106.81415022205296 , _TWO_PI ): m= -1.421e-14 } } else // modulo range: (y..0] { if (m<=y) // Mod(1e-16 , -360. ): m= -360. return 0; if (m>0 ) { if (y+m == y) return 0 ; // just in case... else return y+m; // Mod(-106.81415022205296, -_TWO_PI): m= 1.421e-14 } } return m; } // wrap [rad] angle to [-PI..PI) inline double WrapPosNegPI(double fAng) { return Mod(fAng + _PI, _TWO_PI) - _PI; } // wrap [rad] angle to [0..TWO_PI) inline double WrapTwoPI(double fAng) { return Mod(fAng, _TWO_PI); } // wrap [deg] angle to [-180..180) inline double WrapPosNeg180(double fAng) { return Mod(fAng + 180., 360.) - 180.; } // wrap [deg] angle to [0..360) inline double Wrap360(double fAng) { return Mod(fAng ,360.); } 

Solution à un temps constant one-liner:

D’accord, c’est une double couche si vous comptez la deuxième fonction pour la forme [min,max) , mais suffisamment proche – vous pouvez les fusionner quand même.

 /* change to `float/fmodf` or `long double/fmodl` or `int/%` as appropriate */ /* wrap x -> [0,max) */ double wrapMax(double x, double max) { /* integer math: `(max + x % max) % max` */ return fmod(max + fmod(x, max), max); } /* wrap x -> [min,max) */ double wrapMinMax(double x, double min, double max) { return min + wrapMax(x - min, max - min); } 

Ensuite, vous pouvez simplement utiliser deltaPhase = wrapMinMax(deltaPhase, -M_PI, +M_PI) .

La solution est à temps constant, ce qui signifie que le temps nécessaire ne dépend pas de la distance qui [-PI,+PI) votre valeur de [-PI,+PI) – pour le meilleur ou pour le pire.

Vérification:

Maintenant, je ne m’attends pas à ce que vous me croyiez sur parole, alors voici quelques exemples, y compris les conditions aux limites. J’utilise des entiers pour plus de clarté, mais cela fonctionne à peu près de la même manière avec fmod() et float:

  • Positif x :
    • wrapMax(3, 5) == 3 : (5 + 3 % 5) % 5 == (5 + 3) % 5 == 8 % 5 == 3
    • wrapMax(6, 5) == 1 : (5 + 6 % 5) % 5 == (5 + 1) % 5 == 6 % 5 == 1
  • Négatif x :
    • Remarque: Ceux ci supposent qu’un entier modulo copie le signe de la main gauche; sinon, vous obtenez le cas ci-dessus (“positif”).
    • wrapMax(-3, 5) == 2 : (5 + (-3) % 5) % 5 == (5 - 3) % 5 == 2 % 5 == 2
    • wrapMax(-6, 5) == 4 : (5 + (-6) % 5) % 5 == (5 - 1) % 5 == 4 % 5 == 4
  • Limites:
    • wrapMax(0, 5) == 0 : (5 + 0 % 5) % 5 == (5 + 0) % 5 == 5 % 5 == 0
    • wrapMax(5, 5) == 0 : (5 + 5 % 5) % 5 == (5 + 0) % 5== 5 % 5 == 0
    • wrapMax(-5, 5) == 0 : (5 + (-5) % 5) % 5 == (5 + 0) % 5 == 5 % 5 == 0
      • Remarque: éventuellement -0 au lieu de +0 pour la virgule flottante.

La fonction wrapMinMax fonctionne à peu près de la même façon: insérer x dans [min,max) revient à utiliser x - min dans [0,max-min) , puis (ré) append min au résultat.

Je ne sais pas ce qui se passerait avec un max négatif, mais n’hésitez pas à le vérifier vous-même!

Il existe également fmod fonction fmod dans math.h mais le signe pose des problèmes, de sorte qu’une opération ultérieure est nécessaire pour que le résultat soit dans la plage appropriée (comme vous le faites déjà avec le fmod while). Pour les grandes valeurs de deltaPhase cela est probablement plus rapide que de soustraire / append `M_TWOPI ‘des centaines de fois.

 deltaPhase = fmod(deltaPhase, M_TWOPI); 

EDIT: Je n’ai pas essayé intensément, mais je pense que vous pouvez utiliser fmod cette façon en gérant les valeurs positives et négatives différemment:

  if (deltaPhase>0) deltaPhase = fmod(deltaPhase+M_PI, 2.0*M_PI)-M_PI; else deltaPhase = fmod(deltaPhase-M_PI, 2.0*M_PI)+M_PI; 

Le temps de calcul est constant (contrairement à la solution while qui ralentit à mesure que la valeur absolue de deltaPhase augmente)

Si jamais votre angle d’entrée peut atteindre des valeurs arbitrairement élevées, et si la continuité est importante, vous pouvez également essayer

 atan2(sin(x),cos(x)) 

Cela préservera mieux la continuité de sin (x) et cos (x) que modulo pour les valeurs élevées de x, en particulier en simple précision (float).

En effet, exact_value_of_pi – double_precision_approximation ~ = 1.22e-16

Par ailleurs, la plupart des bibliothèques / matériels utilisent une approximation de haute précision de PI pour appliquer le modulo lors de l’évaluation des fonctions sortinggonomésortingques (bien que la famille x86 en utilise une assez pauvre).

Le résultat peut être dans [-pi, pi], vous devrez vérifier les limites exactes.

Personnellement, j’empêcherais n’importe quel angle d’atteindre plusieurs révolutions en enveloppant systématiquement et en collant à une solution fmod comme celle de boost.

Je ferais ceci:

 double wrap(double x) { return x-2*M_PI*floor(x/(2*M_PI)+0.5); } 

Il y aura des erreurs numériques importantes. La meilleure solution aux erreurs numériques consiste à stocker votre phase mise à l’échelle de 1 / PI ou de 1 / (2 * PI) et, en fonction de vos activités, de les stocker en tant que point fixe.

Au lieu de travailler en radians, utilisez des angles réduits de 1 / (2π) et utilisez modf, floor, etc. Reconvertissez en radians pour utiliser les fonctions de bibliothèque.

Cela a également pour effet que tourner dix mille et demi tours équivaut à faire tourner demi puis dix mille tours, ce qui n’est pas garanti si vos angles sont exprimés en radians, car vous avez une représentation exacte dans la valeur à virgule flottante plutôt qu’une sum approximative. représentations:

 #include  #include  float wrap_rads ( float r ) { while ( r > M_PI ) { r -= 2 * M_PI; } while ( r <= -M_PI ) { r += 2 * M_PI; } return r; } float wrap_grads ( float r ) { float i; r = modff ( r, &i ); if ( r > 0.5 ) r -= 1; if ( r <= -0.5 ) r += 1; return r; } int main () { for (int rotations = 1; rotations < 100000; rotations *= 10 ) { { float pi = ( float ) M_PI; float two_pi = 2 * pi; float a = pi; a += rotations * two_pi; std::cout << rotations << " and a half rotations in radians " << a << " => " << wrap_rads ( a ) / two_pi << '\n' ; } { float pi = ( float ) 0.5; float two_pi = 2 * pi; float a = pi; a += rotations * two_pi; std::cout << rotations << " and a half rotations in grads " << a << " => " << wrap_grads ( a ) / two_pi << '\n' ; } std::cout << '\n'; }} 

J’ai rencontré cette question lorsque j’ai cherché comment insérer une valeur à virgule flottante (ou un double) entre deux nombres arbitraires. Il ne répondait pas spécifiquement à mon cas et j’ai donc élaboré ma propre solution, que l’on peut voir ici. Cela prendra une valeur donnée et l’enroulera entre lowerBound et upperBound, où upperBound se mariera parfaitement à lowerBound, de sorte qu’ils soient équivalents (c’est-à-dire: 360 degrés == 0 degré pour que 360 ​​soit renvoyé à 0).

Espérons que cette réponse est utile pour les autres qui tombent sur cette question à la recherche d’une solution englobante plus générique.

 double boundBetween(double val, double lowerBound, double upperBound){ if(lowerBound > upperBound){std::swap(lowerBound, upperBound);} val-=lowerBound; //adjust to 0 double rangeSize = upperBound - lowerBound; if(rangeSize == 0){return upperBound;} //avoid dividing by 0 return val - (rangeSize * std::floor(val/rangeSize)) + lowerBound; } 

Une question connexe pour les entiers est disponible ici: Algorithme propre et efficace pour le wrapping d’entiers en C ++

Voici une version pour les autres personnes trouvant cette question pouvant utiliser le C ++ avec Boost:

 #include  #include  template inline T normalizeRadiansPiToMinusPi(T rad) { // copy the sign of the value in radians to the value of pi T signedPI = boost::math::copysign(boost::math::constants::pi(),rad); // set the value of rad to the appropriate signed value between pi and -pi rad = fmod(rad+signedPI,(2*boost::math::constants::pi())) - signedPI; return rad; } 

Version C ++ 11, pas de dépendance Boost:

 #include  // Bring the 'difference' between two angles into [-pi; pi]. template  T normalizeRadiansPiToMinusPi(T rad) { // Copy the sign of the value in radians to the value of pi. T signed_pi = std::copysign(M_PI,rad); // Set the value of difference to the appropriate signed value between pi and -pi. rad = std::fmod(rad + signed_pi,(2 * M_PI)) - signed_pi; return rad; } 

Dans le cas où fmod () est implémenté par division tronquée et a le même signe que le dividende , il peut être utilisé pour résoudre le problème général de la manière suivante:

Pour le cas de (-PI, PI]:

 if (x > 0) x = x - 2PI * ceil(x/2PI) #Shift to the negative regime return fmod(x - PI, 2PI) + PI 

Et pour le cas de [-PI, PI):

 if (x < 0) x = x - 2PI * floor(x/2PI) #Shift to the positive regime return fmod(x + PI, 2PI) - PI 

[Notez qu'il s'agit d'un pseudocode; mon original était écrit en Tcl et je ne voulais pas torturer tout le monde avec ça. J'avais besoin du premier cas, donc dû comprendre cela.]

Une solution testée à deux couches, non itérative, pour la normalisation d’angles arbitraires à [-π, π):

 double normalizeAngle(double angle) { double a = fmod(angle + M_PI, 2 * M_PI); return a >= 0 ? (a - M_PI) : (a + M_PI); } 

De même, pour [0, 2π):

 double normalizeAngle(double angle) { double a = fmod(angle, 2 * M_PI); return a >= 0 ? a : (a + 2 * M_PI); } 

deltaPhase -= floor(deltaPhase/M_TWOPI)*M_TWOPI;

Le moyen que vous avez suggéré est le meilleur. C’est le plus rapide pour les petites déviations. Si les angles de votre programme sont constamment déviés dans la plage appropriée, vous ne devriez que rarement rencontrer de grandes valeurs en dehors de la plage. Par conséquent, payer le coût d’un code arithmétique modulaire complexe à chaque tour semble inutile. Les comparaisons sont peu coûteuses comparées à l’arithmétique modulaire ( http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-utiliser-the-modulus-operator-with-caution/ ).

En C99:

 float unwindRadians( float radians ) { const bool radiansNeedUnwinding = radians < -M_PI || M_PI <= radians; if ( radiansNeedUnwinding ) { if ( signbit( radians ) ) { radians = -fmodf( -radians + M_PI, 2.f * M_PI ) + M_PI; } else { radians = fmodf( radians + M_PI, 2.f * M_PI ) - M_PI; } } return radians; } 

Si vous établissez un lien avec la bibliothèque glibc (y compris l’implémentation de newlib), vous pouvez accéder aux fonctions privées __ieee754_rem_pio2f () et __ieee754_rem_pio2 ():

 extern __int32_t __ieee754_rem_pio2f (float,float*); float wrapToPI(float xf){ const float p[4]={0,M_PI_2,M_PI,-M_PI_2}; float yf[2]; int q; int qmod4; q=__ieee754_rem_pio2f(xf,yf); /* xf = q * M_PI_2 + yf[0] + yf[1] / * yf[1] << y[0], not sure if it could be ignored */ qmod4= q % 4; if (qmod4==2) /* (yf[0] > 0) defines interval (-pi,pi]*/ return ( (yf[0] > 0) ? -p[2] : p[2] ) + yf[0] + yf[1]; else return p[qmod4] + yf[0] + yf[1]; } 

Edit: Je viens de me rendre compte que vous devez créer un lien vers libm.a, je ne pouvais pas trouver les symboles déclarés dans libm.so

J’ai utilisé (en python):

 def WrapAngle(Wrapped, UnWrapped ): TWOPI = math.pi * 2 TWOPIINV = 1.0 / TWOPI return UnWrapped + round((Wrapped - UnWrapped) * TWOPIINV) * TWOPI 

équivalent c-code:

 #define TWOPI 6.28318531 double WrapAngle(const double dWrapped, const double dUnWrapped ) { const double TWOPIINV = 1.0/ TWOPI; return dUnWrapped + round((dWrapped - dUnWrapped) * TWOPIINV) * TWOPI; } 

remarquez que cela le place dans le domaine encapsulé de +/- 2pi, donc pour le domaine +/- pi, vous devez le gérer ultérieurement, par exemple:

 if( angle > pi): angle -= 2*math.pi