Я ищу хороший C-код, который будет эффективно реализован:
while (deltaPhase >= M_PI) deltaPhase -= M_TWOPI;
while (deltaPhase < -M_PI) deltaPhase += M_TWOPI;
Каковы мои варианты?
Я ищу хороший C-код, который будет эффективно реализован:
while (deltaPhase >= M_PI) deltaPhase -= M_TWOPI;
while (deltaPhase < -M_PI) deltaPhase += M_TWOPI;
Каковы мои варианты?
Изменить 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 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.);
}
Хорошо, это двухстрочный, если вы считаете вторую функцию для формы [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
работает примерно так: wrapping x
to [min,max)
совпадает с оберткой x - min
до [0,max-min)
, а затем (повторно) добавляет min
к результату.
Я не знаю, что произойдет с отрицательным максимумом, но не стесняйтесь проверить это сами!
Существует также функция fmod
в math.h
, но знак вызывает проблему, так что требуется следующая операция, чтобы сделать результат fir в правильном диапазоне (например, вы уже делаете с while). Для больших значений deltaPhase
это, вероятно, быстрее, чем вычитание/добавление `M_TWOPI 'сотни раз.
deltaPhase = fmod(deltaPhase, M_TWOPI);
EDIT:
Я не пытался интенсивно, но я думаю, что вы можете использовать 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)
Если ваш входной угол может достигать сколь угодно высоких значений, и если непрерывность имеет значение, вы также можете попробовать
atan2(sin(x),cos(x))
Это сохранит непрерывность sin (x) и cos (x) лучше, чем по модулю при высоких значениях x, особенно в одиночной точности (float).
Действительно, exact_value_of_pi - double_precision_approximation ~ = 1.22e-16
С другой стороны, большинство библиотек/аппаратных средств используют высокоточную аппроксимацию PI для применения модуля при оценке тригонометрических функций (хотя известно, что семейство x86 использует довольно бедный).
Результат может быть в [-pi, pi], вам нужно будет проверить точные границы.
Персоналий, я бы предотвратил любой угол, чтобы достичь нескольких оборотов, обертываясь систематически и придерживаясь решения fmod, такого как boost.
Я бы сделал это:
double wrap(double x) {
return x-2*M_PI*floor(x/(2*M_PI)+0.5);
}
Будут значительные числовые ошибки. Лучшим решением для числовых ошибок является сохранение фазы с шагом 1/PI или 1/(2 * PI), и в зависимости от того, что вы делаете, сохраняете их в качестве фиксированной точки.
Вместо работы в радианах используйте углы, масштабируемые 1/(2 & pi;), и используйте modf, пол и т.д. Преобразуйте обратно в радианы, чтобы использовать библиотечные функции.
Это также приводит к тому, что вращение в десять тысяч с половиной оборотов совпадает с вращающейся половиной, а затем с десятью тысячами оборотов, что не гарантируется, если ваши углы находятся в радианах, так как у вас есть точное представление в значении с плавающей запятой чем суммирование приближенных представлений:
#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';
}}
Я столкнулся с этим вопросом при поиске того, как обернуть значение с плавающей запятой (или двойным) между двумя произвольными числами. Он не отвечал конкретно за мой случай, поэтому я разработал свое собственное решение, которое можно увидеть здесь. Это займет определенное значение и обернет его между 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;
}
Соответствующий вопрос для целых чисел доступен здесь: Чистый, эффективный алгоритм для объединения целых чисел в С++
Вот версия для других людей, которые находят этот вопрос, который может использовать С++ с 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;
}
Версия С++ 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;
}
В случае, когда fmod() реализуется через усеченное деление и имеет тот же знак, что и дивиденд, его можно использовать в решить общую проблему:
В случае (-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, и я не хотел мучить всех этим. Мне нужен был первый случай, поэтому пришлось это выяснить.]
Двухстрочное, неитеративное, проверенное решение для нормализации произвольных углов к [-π, π):
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);
}
deltaPhase -= floor(deltaPhase/M_TWOPI)*M_TWOPI;
Я использовал (в python):
def WrapAngle(Wrapped, UnWrapped ):
TWOPI = math.pi * 2
TWOPIINV = 1.0 / TWOPI
return UnWrapped + round((Wrapped - UnWrapped) * TWOPIINV) * TWOPI
эквивалент 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;
}
обратите внимание, что это приведет его к завернутому домену +/- 2pi, так что для домена +/- pi вам необходимо обработать это потом:
if( angle > pi):
angle -= 2*math.pi
Рекомендуемый вам способ лучше всего. Это самый быстрый для небольших прогибов. Если углы в вашей программе постоянно отклоняются в надлежащий диапазон, вы должны редко встречаться в больших значениях вне диапазона. Поэтому оплата стоимости сложного модульного арифметического кода каждый раунд кажется расточительной. Сравнение дешево по сравнению с модульной арифметикой (http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-use-the-modulus-operator-with-caution/).
В 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;
}
Если вы ссылаетесь на glibc libm (включая реализацию newlib), вы можете получить доступ __ieee754_rem_pio2f() и __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];
}
Изменить: только что поняли, что вам нужно установить ссылку на libm.a, я не смог найти символы, объявленные в libm.so