私は効果的に達成するいくつかの素晴らしい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;
私のオプションは何ですか?
2013 年 4 月 19 日編集:
aka.nice と arr_sea で指摘されているように、境界ケースを処理するように Modulo 関数が更新されました。
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.);
}
入力角度が任意に高い値に達する可能性があり、連続性が重要な場合は、試すこともできます
atan2(sin(x),cos(x))
これにより、sin(x) と cos(x) の連続性が、特に単精度 (float) で、x の値が大きい場合に modulo よりも良好に保持されます。
確かに、exact_value_of_pi - double_precision_approximation ~= 1.22e-16
一方、ほとんどのライブラリ/ハードウェアは、三角関数を評価するときにモジュロを適用するために PI の高精度近似を使用します (ただし、x86 ファミリはかなり貧弱なものを使用することが知られています)。
結果は [-pi,pi] にある可能性があります。正確な境界を確認する必要があります。
個人的には、体系的にラップすることで角度が数回転に達するのを防ぎ、ブーストのような fmod ソリューションに固執します。
fmod
関数 inもありますmath.h
が、符号によって問題が発生するため、結果を適切な範囲にするために後続の操作が必要になります (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;
計算時間は一定です (deltaPhase の絶対値が増加するにつれて遅くなる while ソリューションとは異なります)。
私はこれをします:
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 などを使用します。ラジアンに戻してライブラリ関数を使用します。
これには、1 万 5 回転の回転が 1 万回転の 2 分の 1 回転と同じであるという効果もあります。これは、角度がラジアンの場合は保証されません。表現:
#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';
}}
2つの任意の数値の間で浮動小数点値(またはdouble)をラップする方法を検索しているときに、この質問に遭遇しました。私の場合は特に答えられなかったので、ここで見ることができる独自の解決策を考え出しました。これは、指定された値を取り、それを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++で整数をラップするためのクリーンで効率的なアルゴリズム
BoostでC++を使用できる、この質問を見つけた他の人向けのバージョンを次に示します。
#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;
}
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 で書かれていて、それでみんなを苦しめたくありませんでした。最初のケースが必要だったので、これを理解しなければなりませんでした。]
任意の角度を [-π, π) に正規化するための 2 ライナーの非反復テスト済みソリューション:
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/effective-c-tip-13-use-the-modulus-operator-with-caution/ )。