C：如何将浮点数包装到区间 [-pi, pi)
Posted
技术标签:
【中文标题】C：如何将浮点数包装到区间 [-pi, pi)【英文标题】：C: How to wrap a float to the interval [-pi, pi) 【发布时间】：2011-06-05 17:06:51 【问题描述】：我正在寻找一些可以有效完成的不错的 C 代码：
while (deltaPhase >= M_PI) deltaPhase -= M_TWOPI;
while (deltaPhase < -M_PI) deltaPhase += M_TWOPI;
我有什么选择？
【问题讨论】：
不要忘记，while 不仅会累积错误，而且如果输入非常高的值（如 ldexp(M_PI,55)）可能会变成无限循环 【参考方案1】：deltaPhase -= floor(deltaPhase/M_TWOPI)*M_TWOPI;
【讨论】：
此代码会产生不准确的结果，不应使用。fmod 存在是有原因的。
出于好奇，这有什么问题？我已经对其进行了测试，并且效果很好。有人可以举一个失败的例子吗？
鉴于没有人愿意指出其中的缺陷，我投票将其设为 0
@P i: 让 c 在[0, 1) 中，让deltaPhase=-c*PI。然后我们得到-c*PI - (-1)*2*PI，它等于(2-c)*PI，它不在[-pi, pi)中。所以我会撤回你的赞成票。【参考方案2】：
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 的绝对值增加而变慢）
【讨论】：
你能否给出处理负值的完整实现？ 我忘了说如果你使用 gcc，你必须链接到libm.a。
负数不是在 (-pi,pi] 吗？
简单有效，正常使用。【参考方案3】：
不要使用弧度，而是使用按 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';
【讨论】：
太糟糕了，没有标准库函数假设一个完整的圆代表值 1、2 或 4 [取决于是否要计算旋转、圆周率或象限]，因为范围缩小会容易得多。从实际的角度来看，即使在计算之前乘以 pi 会增加一个潜在的舍入误差，但使用每转 2 的幂几乎总是会提高早期计算的精度，而不是最小的舍入误差会损害它。【参考方案4】：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<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.);
【讨论】：
试试这个，在 IEEE 754 双精度中可能会失败（没有扩展精度提升，-ffloat-store）assert(WrapPosNegPI(103.67255756846316) >= - _PI);我找到了以下 Smalltalk sn-p (1 to: 11111 by: 2) 的示例检测：[:i | ((i *Float pi) 前身 / Float pi) floor = i] 一个问题：Mod(x,360.0) 应该将事物包装在 [0,360) 范围内。但是，当所需的返回值为 0.0 时，此 Mod(-1e-16, 360.0) 的实现返回 360.0。这是因为数学试图返回 359.9999999999999999 但不能用双精度表示，因此四舍五入为 360.0。一种解决方法可能是首先插入行“x += 10.0*y;”在 Mod 函数的开头，以防止精度损失导致此故障。肮脏或优雅......你决定:) -1。 方式过于复杂，有很多分支，使用保留标识符（以_[A-Z] 开头的标识符），但也许更重要的是 --- 问题被标记为 C，答案是 C++。
这个Mod()有没有比标准fmod()更好的特定方式？
@Dolda2000：没有好坏之分，只是定义不同而已。有关浮点模函数的可能定义，请参见第 4 节 here。【参考方案5】：
我在搜索如何在两个任意数字之间包装浮点值（或双精度值）时遇到了这个问题。它没有专门针对我的情况回答，所以我制定了自己的解决方案，可以在这里看到。这将取一个给定的值并将其包裹在 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;
此处提供了一个与整数相关的问题： Clean, efficient algorithm for wrapping integers in C++
【讨论】：
【参考方案6】：我用过（在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
【讨论】：
等一下，TWOPI = math.pi？我猜你错过了*2。【参考方案7】：如果 fmod() 是通过截除实现的，并且与dividend 具有相同的符号，则可以利用它来解决一般问题：
对于(-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 编写的，我不想用它来折磨每个人。我需要第一个案例，所以必须弄清楚这一点。]
【讨论】：
@Pavel Ognev 不要对人们的答案做出根本性的改变。【参考方案8】：这是其他人发现可以使用 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;
【讨论】：
在 (-pi/2, pi/2) 之间换一个角度怎么样？ @CroCo 只是把这个函数的输出除以 2，对吧？ 虽然上述方法很简洁，但我只想指出输出角度在 [-pi, pi] 范围内，而不是原始问题所要求的 [-p, pi) 范围内。 哇，我多年来一直在使用这两个版本，但我从未注意到这一点，感谢您的洞察力！在我的用例中这不是问题，我认为针对该确切值的额外 if 语句可以解决它，但我愿意接受更好的方法。boost::math::constants::pi<T>()我的上帝，boost应该真的死了。你必须有一种特殊的才能，使简单的东西在阅读时难以记忆、使用和理解。我知道这是做事的“C++ 方式”，但这意味着 C++ 在此过程中出现了问题。我很高兴我一直避免使用 boost。【参考方案9】：
如果您的输入角度可以达到任意高值，并且如果连续性很重要，您也可以尝试
atan2(sin(x),cos(x))
这将比模数更好地保持 sin(x) 和 cos(x) 的连续性，对于 x 的高值，尤其是在单精度（浮点数）中。
确实，exact_value_of_pi - double_precision_approximation ~= 1.22e-16
另一方面，大多数库/硬件在评估三角函数时使用 PI 的高精度近似值来应用模数（尽管已知 x86 系列使用了一个相当差的函数）。
结果可能在 [-pi,pi] 中，您必须检查确切的界限。
就个人而言，我会通过系统地环绕并坚持使用像 boost 那样的 fmod 解决方案来防止任何角度达到几圈。
【讨论】：
一个聪明的主意，即使你最终没有得到这个实现，它也是一个测试你自己的好方法。太好了！【参考方案10】：我会这样做：
double wrap(double x)
return x-2*M_PI*floor(x/(2*M_PI)+0.5);
会有明显的数字错误。数值误差的最佳解决方案是存储按 1/PI 或 1/(2*PI) 缩放的相位，并根据您的操作将它们存储为固定点。
【讨论】：
【参考方案11】：您建议的方式是最好的。对于小挠度，它是最快的。如果您的程序中的角度不断偏转到适当的范围内，那么您应该很少遇到大的超出范围的值。因此，每轮都支付一个复杂的模算术代码的成本似乎是一种浪费。与模运算 (http://embeddedgurus.com/stack-overflow/2011/02/efficient-c-tip-13-use-the-modulus-operator-with-caution/) 相比，比较便宜。
【讨论】：
【参考方案12】：在 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;
【讨论】：
【参考方案13】：单线恒定时间解：
好的，如果您计算 [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 添加到结果中。
我不知道负最大值会发生什么，但请自行检查！
【讨论】：
真的，你认为fmod 是常数时间吗？ % 甚至不是整数的常数时间。
但这不是我来这里要说的：我要编写一个程序，将您的函数wrapMinMax 应用于随机值。在我写之前，你想打赌wrapMinMax 返回的值低于min 和/或高于max？
@PascalCuoq：嗯，好的，但是执行的指令数量是恒定的。我已经编写了测试这个的程序：codepad.org/aMbhviku --- 请注意它是如何在 10M 测试大小时没有输出的。也就是说，假设max > min。从长远来看，失去一些傲慢/光顾可能会对你有所帮助。
@PascalCuoq：这是整数版本：codepad.org/YQBtMpQZ --- 是的，#define double int 是一个丑陋的黑客，但我很懒。我也忘了把%f改成%d，但是够接近了。
对不起我以前的cmets语气。你的回答让我烦恼的是，我应该简单地表达而不是含糊其辞的讽刺，这是一个非常古老的问题，有很多答案，我真的看不出你的答案如何改进现有的答案，而且还有存在一个正确的答案，但它是一本书的章节，对于这个老问题，目前总结起来工作量太大。本书章节的标题是“三角函数的参数缩减”，其中包含“Payne”和“Hanek”这两个词。【参考方案14】：
如果链接到 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 中声明的符号
【讨论】：
【参考方案15】：将任意角度归一化为 [-π, π) 的双线性非迭代测试解决方案：
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);
【讨论】：
以上是关于C：如何将浮点数包装到区间 [-pi, pi)的主要内容，如果未能解决你的问题，请参考以下文章