C: Cómo envolver un flotador en el intervalo [-pi, pi)

Estoy buscando un buen código C que logre de manera efectiva:

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

¿Cuáles son mis opciones?

Editar 19 de abril de 2013:

La función de módulo se actualizó para manejar casos de límite como notaron aka.nice y 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.); } 

Solución de una línea de tiempo constante:

De acuerdo, es un juego de dos líneas si cuentas la segunda función para la forma [min,max) , pero lo suficientemente cerca, podrías unirlas de todas formas.

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

Entonces, simplemente puede usar deltaPhase = wrapMinMax(deltaPhase, -M_PI, +M_PI) .

Las soluciones son de tiempo constante, lo que significa que el tiempo que toma no depende de cuán lejos esté su valor de [-PI,+PI) – para bien o para mal.

Verificación:

Ahora, no espero que asums mi palabra, así que aquí hay algunos ejemplos, incluidas las condiciones de contorno. Estoy usando números enteros para mayor claridad, pero funciona muy fmod() con fmod() y floats:

  • Positivo 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
  • Negativo x :
    • Nota: Suponen que el módulo entero copia el signo de la mano izquierda; si no, obtienes el caso anterior (“Positivo”).
    • 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
  • Límites:
    • 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
      • Nota: Posiblemente -0 lugar de +0 para coma flotante.

La función wrapMinMax funciona de manera muy parecida: ajustar x a [min,max) es lo mismo que ajustar x - min a [0,max-min) y luego (re) agregar min al resultado.

No sé qué pasaría con un máximo negativo, ¡pero siéntete libre de comprobarlo tú mismo!

También existe la función fmod en math.h pero el signo causa problemas, por lo que se necesita una operación posterior para que el resultado se encuentre en el rango correcto (como lo hace con los while). Para valores grandes de deltaPhase esto es probablemente más rápido que restar / agregar `M_TWOPI ‘cientos de veces.

 deltaPhase = fmod(deltaPhase, M_TWOPI); 

EDITAR: No lo intenté intensamente, pero creo que puedes usar fmod esta forma al manejar los valores positivos y negativos de manera diferente:

  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; 

El tiempo computacional es constante (a diferencia de la solución while que se vuelve más lenta a medida que aumenta el valor absoluto de deltaPhase)

Si alguna vez tu ángulo de entrada puede alcanzar valores arbitrariamente altos, y si la continuidad importa, también puedes probar

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

Esto preservará la continuidad de sin (x) y cos (x) mejor que módulo para valores altos de x, especialmente en precisión simple (float).

De hecho, exact_value_of_pi – double_precision_approximation ~ = 1.22e-16

Por otro lado, la mayoría de las bibliotecas / hardware usan una aproximación de alta precisión de PI para aplicar el módulo cuando se evalúan funciones trigonométricas (aunque se sabe que la familia x86 usa una bastante pobre).

El resultado puede estar en [-pi, pi], tendrás que verificar los límites exactos.

Personalmente, evitaría que cualquier ángulo alcanzara varias revoluciones envolviéndolo sistemáticamente y aferrándome a una solución fmod como la de boost.

Yo haría esto:

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

Habrá errores numéricos significativos. La mejor solución para los errores numéricos es almacenar su fase escalada en 1 / PI o en 1 / (2 * PI) y, dependiendo de lo que esté haciendo, guárdelos como punto fijo.

En lugar de trabajar en radianes, use angularjs escalados en 1 / (2π) y use modf, piso, etc. Convierta de nuevo a radianes para usar las funciones de la biblioteca.

Esto también tiene el efecto de que girar diez mil revoluciones es lo mismo que rotar la mitad que diez mil revoluciones, lo que no está garantizado si tus angularjs están en radianes, ya que tienes una representación exacta en el valor del punto flotante en lugar de sumr aproximadamente representaciones:

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

Encontré esta pregunta cuando buscaba cómo ajustar un valor de coma flotante (o un doble) entre dos números arbitrarios. No respondía específicamente para mi caso, así que resolví mi propia solución que se puede ver aquí. Esto tomará un valor dado y lo envolverá entre lowerBound y upperBound donde upperBound se encuentra perfectamente con lowerBound de modo que sean equivalentes (es decir: 360 grados == 0 grados para que 360 ​​se ajuste a 0)

Esperemos que esta respuesta sea útil para que otros encuentren esta pregunta en busca de una solución límite más genérica.

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

Una pregunta relacionada para enteros está disponible aquí: Algoritmo limpio y eficiente para envolver enteros en C ++

Aquí hay una versión para otras personas que encuentran esta pregunta que puede usar C ++ con 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; } 

Versión de C ++ 11, sin dependencia de 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; } 

En el caso donde fmod () se implementa a través de una división truncada y tiene el mismo signo que el dividendo , se puede aprovechar para resolver el problema general de esta manera:

Para el caso de (-PI, PI):

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

Y para el caso de [-PI, PI):

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

[Tenga en cuenta que esto es pseudocódigo; mi original fue escrito en Tcl, y no quería torturar a todos con eso. Necesitaba el primer caso, así que tuve que resolverlo.]

Una solución de dos líneas, no iterativa, probada para normalizar angularjs arbitrarios a [-π, π):

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

Del mismo modo, para [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;

La forma sugerida que sugeriste es la mejor. Es más rápido para pequeñas deflexiones. Si los angularjs de su progtwig se desvían constantemente dentro del rango adecuado, entonces solo se encontrarán en grandes valores fuera de rango raramente. Por lo tanto, pagar el costo de un complicado código aritmético modular en cada ronda parece un desperdicio. Las comparaciones son baratas en comparación con la aritmética modular ( http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-use-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 enlaza con la libm de glibc (incluida la implementación de newlib), puede acceder a las funciones privadas __ieee754_rem_pio2f () y __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]; } 

Editar: me acabo de dar cuenta de que necesitas vincular a libm.a, no pude encontrar los símbolos declarados en libm.so

He usado (en python):

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

equivalente en código 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; } 

Observe que esto lo trae en el dominio envuelto +/- 2pi, por lo que para el dominio +/- pi debe manejarlo luego como:

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