C: comment envelopper un flotteur dans l'intervalle [- pi, pi)

je suis à la recherche d'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?

28
demandé sur mskfisher 2011-01-08 11:54:54

15 réponses

Edit 19 Avr. 2013:

fonction Modulo mise à jour pour traiter les cas limites tel que noté 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<typename T>
T Mod(T x, T y)
{
    static_assert(!std::numeric_limits<T>::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.);
}
23
répondu Lior Kogan 2013-04-19 20:33:44

One-liner constante de temps de la solution:

D'accord, c'est un double-liner si vous comptez la deuxième fonction pour [min,max) forme, mais assez proche - vous pourriez les fusionner de toute façon.

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

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

les solutions sont temps constant, ce qui signifie que le temps qu'il prend ne dépend pas de la distance de votre valeur est de [-PI,+PI) -pour le meilleur ou pour le pire.

"1519290920 la" Vérification":

maintenant, je ne m'attends pas à ce que vous me croyez sur parole, donc voici quelques exemples, y compris les conditions limites. J'utilise des entiers pour la clarté, mais ça fonctionne à peu près de la même façon avec fmod() et flotte:

  • 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 :
    • Note: ceux-ci supposent que le modulo entier copie le signe de gauche; si non, 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
      • Note: peut-être -0 au lieu de +0 pour la virgule flottante.

la fonction wrapMinMax fonctionne à peu près de la même façon: x à [min,max) est le même que l'habillage x - min à [0,max-min) , et puis (re-)ajout de min pour le résultat.

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

10
répondu Tim Čas 2015-04-25 22:04:03

il y a aussi la fonction fmod dans math.h mais le signe cause des problèmes de sorte qu'une opération ultérieure est nécessaire pour faire le premier résultat dans la bonne gamme (comme vous le faites déjà avec le while). Pour les grandes valeurs de deltaPhase c'est probablement plus rapide que de soustraire/ajouter `M_TWOPI' des centaines de fois.

deltaPhase = fmod(deltaPhase, M_TWOPI);

EDIT: Je ne l'ai pas essayé intensivement mais je pense que vous pouvez utiliser fmod de cette façon en manipulant valeurs positives et négatives différentes:

    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 devient plus lente que la valeur absolue de deltaPhase augmente)

8
répondu jdehaan 2011-01-08 19:37:37

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

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

cela préservera la continuité du sin(x) et du cos(x) mieux que modulo pour des valeurs élevées de x, surtout en précision simple (float).

en effet, exact_value_of_pi-double_precision_approximation ~= 1.22 e-16

d'autre part, 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 trigonométriques (bien que la famille x86 est connue pour utiliser un plutôt pauvre).

résultat pourrait être en [- pi, pi], vous aurez à vérifier les limites exactes.

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

7
répondu aka.nice 2012-06-25 06:58:27

je ferais ceci:

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

il y aura d'importantes erreurs numériques. La meilleure solution aux erreurs numériques est de stocker votre phase graduée par 1/PI ou par 1/(2*PI) et en fonction de ce que vous faites les stocker comme point fixe.

7
répondu Craig McGillivary 2012-10-19 22:06:32

au lieu de travailler en radians, utilisez des angles gradués par 1/(2π) et utilisez modf, sol, etc. Convertissez de nouveau en radians pour utiliser les fonctions de la bibliothèque.

cela a aussi l'effet que la rotation de dix mille et demi révolutions est la même que la rotation de la moitié de dix mille révolutions, ce qui n'est pas garanti si vos angles sont en radians, car vous avez une représentation exacte dans la valeur de la pointe flottante plutôt que de Sommer approximatif représentations:

#include <iostream>
#include <cmath>

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';
}}
6
répondu Pete Kirkham 2011-01-08 11:42:08

j'ai rencontré cette question en cherchant comment envelopper une valeur à virgule flottante (ou un double) entre deux nombres arbitraires. Il n'a pas répondu spécifiquement pour mon cas, donc j'ai travaillé sur ma propre solution qui peut être vu ici. Cela prendra une valeur donnée et l'envelopper entre lowerBound et upperBound où upperBound rencontre parfaitement lowerBound tels qu'ils sont équivalents (c'est-à-dire: 360 degrés == 0 degrés donc 360 s'envelopperait à 0)

j'espère que cette réponse est utile aux autres qui se heurtent à cette question à la recherche d'une solution de délimitation 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 envelopper les entiers en C++

4
répondu M2tM 2017-05-23 11:47:11

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

#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/sign.hpp>

template<typename T>
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<T>(),rad);
  // set the value of rad to the appropriate signed value between pi and -pi
  rad = fmod(rad+signedPI,(2*boost::math::constants::pi<T>())) - signedPI;

  return rad;
} 

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

#include <cmath>

// Bring the 'difference' between two angles into [-pi; pi].
template <typename T>
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;
}
4
répondu Andrew Hundt 2017-01-31 03:12:19

dans le cas où fmod() est mis en œuvre 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 ainsi:

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 que c'est un pseudocode; mon original a été écrit en Tcl, et je ne voulais pas torturer tout le monde avec ça. Je besoin le premier cas, j'ai donc dû le comprendre.]

1
répondu Glenn 2012-04-17 23:56:47

Un deux-liner, non itératif, solution testée pour la normalisation arbitraire des angles de [-π, π):

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);
}
1
répondu mojuba 2016-09-28 16:17:31

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

0
répondu lijie 2011-01-08 09:04:14

la façon suggérée est la meilleure. Il est le plus rapide pour les petites déflexions. Si les angles de votre programme sont constamment déviés dans la bonne plage, alors vous ne devriez rencontrer de grandes valeurs hors de la plage que rarement. Par conséquent, payer le coût d'un code arithmétique modulaire complexe chaque tour semble un gaspillage. Les comparaisons sont peu coûteuses par rapport à l'arithmétique modulaire ( http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-use-the-modulus-operator-with-caution / ).

0
répondu Tom Larkworthy 2013-04-16 11:11:51

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;
}
0
répondu otto 2013-06-25 15:00:05

si vous vous liez avec la libm de glibc (y compris l'implémentation de newlib), vous pouvez accéder __ieee754_rem_pio2f() et __ _ _ ieee754_rem_pio2() les fonctions privées:

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: Viens de réaliser que vous avez besoin de lien pour libm.a, je n'ai pas pu trouver les symboles déclarés dans libm.so

0
répondu xvan 2015-10-03 17:05:46

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-code c:

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

notez que cela le ramène dans le domaine enveloppé + / - 2pi donc pour le domaine + / - pi vous devez gérer cela après comme:

if( angle > pi):
    angle -= 2*math.pi
0
répondu Henrik 2015-11-17 14:19:38