C: Как обернуть поплавок в интервал [-pi, pi)

Question

Я ищу хороший C-код, который будет эффективно выполнять:

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

Какие у меня варианты?

Answer 1

Изменить 19 апреля 2013 г .:

Обновлена ​​функция Modulo для обработки граничных случаев, как отмечено aka.nice и 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.);
}

Answer 2

Однострочное решение с постоянным временем:

Ладно, если вы рассчитываете вторую функцию для формы [min,max), это две строки, но достаточно близко - вы можете объединить их вместе в любом случае.

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

Тогда вы можете просто использовать deltaPhase = wrapMinMax(deltaPhase, -M_PI, +M_PI).

Решения являются постоянными, то есть время, которое требуется, не зависит от того, как далеко ваше значение от [-PI,+PI) - в лучшую сторонуили, что еще хуже.

Проверка:

Теперь я не ожидаю, что вы поверите мне на слово, поэтому приведу несколько примеров, включая граничные условия.Я использую целые числа для ясности, но он работает почти так же с fmod() и с плавающей точкой:

• Положительный 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
• Отрицательное x:
  • Примечание: Предполагается, что целое число по модулю копирует знак левой руки;если нет, вы получите вышеуказанный («положительный») случай.
  • 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
• Границы:
  • 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
    • Примечание: Возможно -0 вместо +0 для чисел с плавающей запятой.

Функция wrapMinMax работает практически так же: перенос x до [min,max) аналогичен переносу x - min до [0,max-min), а затем (повторному) добавлению min крезультат.

Я не знаю, что произойдет с отрицательным максимумом, но не стесняйтесь проверить это сами!

Answer 3

В math.h также есть функция fmod, но знак вызывает проблемы, так что необходима дальнейшая операция, чтобы получить результат в нужном диапазоне (как вы уже делаете с while).Для больших значений deltaPhase это, вероятно, быстрее, чем вычитание / добавление `M_TWOPI 'сотни раз.

deltaPhase = fmod(deltaPhase, M_TWOPI);

РЕДАКТИРОВАТЬ: Я не пробовал интенсивно, но я думаю, что выможно использовать fmod таким образом, обрабатывая положительные и отрицательные значения по-разному:

    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;

Время вычислений является постоянным (в отличие от решения while, которое замедляется при увеличении абсолютного значения deltaPhase)

Answer 4

Если когда-либо ваш входной угол может достигать сколь угодно высоких значений, и если непрерывность имеет значение, вы также можете попробовать

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

Это сохранит непрерывность sin (x) и cos (x) лучше, чем по модулю для больших значений x, особенно с одинарной точностью (float).

Действительно, точное_значение_д_двойного_двойного_приложения ~ = 1,22е-16

С другой стороны, большинство библиотек / аппаратных средств используют высокоточное приближение PI для применения модуля по оценке тригонометрических функций (хотя известно, что семейство x86 использует довольно скудную).

Результат может быть в [-pi, pi], вам нужно будет проверить точные границы.

Лично я бы помешал любому углу достичь нескольких оборотов, систематически оборачиваясь и придерживаясь решения fmod, такого как Boost.

Answer 5

Я бы сделал это:

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

Там будут значительные числовые ошибки. Лучшее решение для числовых ошибок - сохранить вашу фазу в масштабе 1 / PI или 1 / (2 * PI) и в зависимости от того, что вы делаете, сохранить их как фиксированную точку.

Answer 6

Вместо работы в радианах, используйте углы, масштабированные на 1 / (2π) и используйте modf, floor и т. Д. Преобразуйте обратно в радианы, чтобы использовать библиотечные функции.

Это также приводит к тому, что вращение на десять тысяч с половиной оборотов такое же, как вращение на половину, а не на десять тысяч оборотов, что не гарантируется, если ваши углы указаны в радианах, поскольку у вас есть точное представление в плавающемзначение точки, а не суммирование приблизительных представлений:

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

Answer 7

Вот версия для других людей, нашедших этот вопрос, которые могут использовать C ++ с 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 версия, без зависимости Boost:

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

Answer 8

Я столкнулся с этим вопросом, когда искал, как обернуть значение с плавающей запятой (или двойное число) между двумя произвольными числами.Он не отвечал конкретно для моего случая, поэтому я разработал собственное решение, которое можно увидеть здесь.Это примет заданное значение и обернет его между lowerBound и upperBound, где upperBound отлично встречается с lowerBound, так что они эквивалентны (то есть: 360 градусов == 0 градусов, поэтому 360 будет переноситься в 0)

Надеюсь, этот ответ полезендругим, кто сталкивается с этим вопросом в поисках более общего ограничивающего решения.

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

Здесь можно найти связанный вопрос для целых чисел: Чистый, эффективный алгоритм для упаковки целых чисел в C ++

Answer 9

Двухлинейное, не итеративное, проверенное решение для нормализации произвольных углов к [-π, π):

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

Аналогично для [0, 2π):

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

Answer 10

В случае, если fmod () реализован через усеченное деление и имеет тот же знак, что и divind , можно воспользоваться преимуществом для решения общей проблемы следующим образом:

Для случая (-PI, PI):

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

А для случая [-PI, PI):

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

[Обратите внимание, что это псевдокод; мой оригинал был написан на Tcl, и я не хотел мучить всех этим. Мне нужен был первый случай, поэтому я должен был это выяснить.]