C:如何将浮动包装到区间

我正在寻找一些有效的C代码:

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

我有什么选择?

编辑2013年4月19日:

更新模数函数以处理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 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.); } 

单线恒定时间解决方案:

好吧,如果你计算[min,max)forms的第二个函数,但它足够接近它是一个双线程 – 你可以将它们合并在一起。

 /* 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)距离 – 无论好坏。

validation:

现在,我不指望你接受我的话,所以这里有一些例子,包括边界条件。 我使用整数来表示清晰度,但它与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添加到结果中。

我不知道最大的负面会发生什么,但请随意检查一下!

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的绝对值增加而变慢)

如果您的输入角度可以达到任意高的值,并且连续性很重要,您也可以尝试

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

对于高x值,这将保持sin(x)和cos(x)的连续性优于模数,特别是在单精度(浮点数)中。

实际上,exact_value_of_pi – double_precision_approximation~ = 1.22e-16

另一方面,大多数库/硬件在评估三角函数时使用PI的高精度近似来应用模(尽管已知x86族使用相当差的函数)。

结果可能在[-pi,pi]中,您必须检查确切的界限。

Personaly,我会阻止任何角度通过系统包装达到几个转数并坚持像一个提升的fmod解决方案。

我会这样做:

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

会有重大的数字错误。 数值误差的最佳解决方案是将相位按1 / PI或1 /(2 * PI)存储,并根据您的操作将它们存储为固定点。

使用以1 /(2π)缩放的角度并使用modf,floor等,而不是以弧度工作。转换回弧度以使用库函数。

这也有旋转一万半转数与旋转一半转一万转相同的效果,如果您的角度是弧度,则无法保证,因为您在浮点值中具有精确表示而不是求和近似值陈述:

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

我在搜索如何在两个任意数字之间包装浮点值(或双精度)时遇到了这个问题。 它没有专门回答我的情况,所以我制定了自己的解决方案,可以在这里看到。 这将采用给定的值并将其包装在lowerBound和upperBound之间,其中up​​perBound完全符合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 ++中包装整数的干净,高效的算法

以下是其他人发现此问题的版本,可以将C ++与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; } 

C ++ 11版本,没有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; } 

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

你建议的方式最好。 对于小偏转来说,它是最快的。 如果程序中的角度不断被偏转到适当的范围内,那么你应该很少遇到大的超出范围的值。 因此,每轮支付复杂的模块化算术代码的成本似乎是浪费的。 与模运算相比,比较便宜( 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中声明的符号

我用过(在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