Commit 1d56057c authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files

Change to power function in FMath, and corresponding changes to FOriginalKernel. Tested

parent 33701cb2
......@@ -119,19 +119,19 @@ class FRotationOriginalKernel : public FAbstractKernels<CellClass,ContainerClass
if( l >= 0 && -l <= k && k <= l && FMath::Abs(k) <= m && m <= l ){
FReal sum = 0;
for(int n = FMath::Max(-(m+k),0) ; n <= FMath::Min(l-m,l-k) ; ++n){
sum += FMath::pow(FReal(-1.0), l-m-n) * combin(l-k, n) * combin(l+k, l-m-n) * FMath::pow(FReal(1.0)+cosTheta,n) * FMath::pow(FReal(1.0)-cosTheta, l-m-n);
sum += FMath::pow(FReal(-1.0), FMath::Abs(l-m-n)) * combin(l-k, n) * combin(l+k, l-m-n) * FMath::pow(FReal(1.0)+cosTheta,FMath::Abs(n)) * FMath::pow(FReal(1.0)-cosTheta, FMath::Abs(l-m-n));
}
return (FReal(1.0)/FMath::pow(FReal(2.0),l)) * FMath::Sqrt((fact(l-m)*fact(l+m))/(fact(l-k)*fact(l+k)))
* FMath::pow(FReal(1.0) + FReal(Sgn(k))*cosTheta, FMath::Abs(k)) * FMath::pow(sinTheta, m - FMath::Abs(k)) * sum;
return (FReal(1.0)/FMath::pow(FReal(2.0),FMath::Abs(l))) * FMath::Sqrt((fact(l-m)*fact(l+m))/(fact(l-k)*fact(l+k)))
* FMath::pow(FReal(1.0) + FReal(Sgn(k))*cosTheta, FMath::Abs(k)) * FMath::pow(sinTheta, FMath::Abs(m - FMath::Abs(k))) * sum;
}
// P^{l}_{m,k} = (-1)^{m+k} P^{l}_{k,m} , l \ge 0, -l \leq m \le 0, |m| \leq k \leq l
else if( (l > 0 && -l <= m && m < 0 && FMath::Abs(m) <= k && k <= l)
|| (l > 0 && 0 <= m && m < l && m < k && k <= l )) {
return FMath::pow(FReal(-1.0), m+k) * d_lmk_analytical(cosTheta,sinTheta, l, k ,m);
return FMath::pow(FReal(-1.0), FMath::Abs(m+k)) * d_lmk_analytical(cosTheta,sinTheta, l, k ,m);
}
// P^{l}_{m,k} = (-1)^{m+k} P^{l}_{k,m} , l \ge 0, 0 \leq m \le l, m \le k \leq l
else if( l > 0 && -l <= m && m < l && -l <= k && k < -m ){
return FMath::pow(FReal(-1.0), m+k) * d_lmk_analytical(cosTheta,sinTheta, l, -m, -k);
return FMath::pow(FReal(-1.0), FMath::Abs(m+k)) * d_lmk_analytical(cosTheta,sinTheta, l, -m, -k);
}
else{
printf("Error that should not be possible!\n");
......
......@@ -21,7 +21,7 @@
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class
* @class This class is a basic implementation of Complexe.
* Please read the license
*
* Propose basic complexe class.
......
......@@ -23,228 +23,234 @@
#include "FGlobal.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class
* Please read the license
*
* Propose basic math functions or indirections to std math.
*/
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class
* Please read the license
*
* Propose basic math functions or indirections to std math.
*/
struct FMath{
static const FReal FPi; //< Pi constant
static const FReal FTwoPi; //< 2 Pi constant
static const FReal FPiDiv2; //< Pi/2 constant
static const FReal Epsilon; //< Epsilon
/** To get absolute value */
template <class NumType>
static NumType Abs(const NumType inV){
return (inV < 0 ? -inV : inV);
}
static const FReal FPi; //< Pi constant
static const FReal FTwoPi; //< 2 Pi constant
static const FReal FPiDiv2; //< Pi/2 constant
static const FReal Epsilon; //< Epsilon
/** To get max between 2 values */
template <class NumType>
static NumType Max(const NumType inV1, const NumType inV2){
return (inV1 > inV2 ? inV1 : inV2);
}
/** To get absolute value */
template <class NumType>
static NumType Abs(const NumType inV){
return (inV < 0 ? -inV : inV);
}
/** To get min between 2 values */
template <class NumType>
static NumType Min(const NumType inV1, const NumType inV2){
return (inV1 < inV2 ? inV1 : inV2);
}
/** To get max between 2 values */
template <class NumType>
static NumType Max(const NumType inV1, const NumType inV2){
return (inV1 > inV2 ? inV1 : inV2);
}
/** To know if 2 values seems to be equal */
template <class NumType>
static bool LookEqual(const NumType inV1, const NumType inV2){
return (Abs(inV1-inV2) < std::numeric_limits<NumType>::epsilon());
//const FReal relTol = FReal(0.00001);
//const FReal absTol = FReal(0.00001);
//return (Abs(inV1 - inV2) <= Max(absTol, relTol * Max(Abs(inV1), Abs(inV2))));
}
/** To get min between 2 values */
template <class NumType>
static NumType Min(const NumType inV1, const NumType inV2){
return (inV1 < inV2 ? inV1 : inV2);
}
/** To know if 2 values seems to be equal */
template <class NumType>
static FReal RelatifDiff(const NumType inV1, const NumType inV2){
return Abs(inV1 - inV2)*Abs(inV1 - inV2)/Max(Abs(inV1*inV1), Abs(inV2*inV2));
}
/** To know if 2 values seems to be equal */
template <class NumType>
static bool LookEqual(const NumType inV1, const NumType inV2){
return (Abs(inV1-inV2) < std::numeric_limits<NumType>::epsilon());
//const FReal relTol = FReal(0.00001);
//const FReal absTol = FReal(0.00001);
//return (Abs(inV1 - inV2) <= Max(absTol, relTol * Max(Abs(inV1), Abs(inV2))));
}
/** To get floor of a FReal */
static float dfloor(const float inValue){
return floorf(inValue);
}
static double dfloor(const double inValue){
return floor(inValue);
}
/** To know if 2 values seems to be equal */
template <class NumType>
static FReal RelatifDiff(const NumType inV1, const NumType inV2){
return Abs(inV1 - inV2)*Abs(inV1 - inV2)/Max(Abs(inV1*inV1), Abs(inV2*inV2));
}
/** To get ceil of a FReal */
static float Ceil(const float inValue){
return ceilf(inValue);
}
static double Ceil(const double inValue){
return ceil(inValue);
}
/** To get floor of a FReal */
static float dfloor(const float inValue){
return floorf(inValue);
}
static double dfloor(const double inValue){
return floor(inValue);
}
/** To get pow */
static double pow(double x, double y){
return ::pow(x,y);
}
static double pow(float x, float y){
return ::powf(x,y);
}
template <class NumType>
static NumType pow(const NumType inValue, int power){
if(power<0) power = -power;
NumType result = 1;
while(power-- > 0) result *= inValue;
return result;
}
/** To get ceil of a FReal */
static float Ceil(const float inValue){
return ceilf(inValue);
}
static double Ceil(const double inValue){
return ceil(inValue);
}
/** To get pow of 2 */
static int pow2(const int power){
return (1 << power);
}
/** To get pow */
static double pow(double x, double y){
return ::pow(x,y);
}
static double pow(float x, float y){
return ::powf(x,y);
}
template <class NumType>
static NumType pow(const NumType inValue, int power){
NumType result = 1;
while(power-- > 0) result *= inValue;
return result;
}
/** To know if a value is between two others */
template <class NumType>
static bool Between(const NumType inValue, const NumType inMin, const NumType inMax){
return ( inMin <= inValue && inValue < inMax );
}
/** To get pow of 2 */
static int pow2(const int power){
return (1 << power);
}
/** To get sqrt of a FReal */
static float Sqrt(const float inValue){
return sqrtf(inValue);
}
static double Sqrt(const double inValue){
return sqrt(inValue);
}
/** To know if a value is between two others */
template <class NumType>
static bool Between(const NumType inValue, const NumType inMin, const NumType inMax){
return ( inMin <= inValue && inValue < inMax );
}
/** To get Log of a FReal */
static float Log(const float inValue){
return logf(inValue);
}
static double Log(const double inValue){
return log(inValue);
}
/** To get sqrt of a FReal */
static float Sqrt(const float inValue){
return sqrtf(inValue);
}
static double Sqrt(const double inValue){
return sqrt(inValue);
}
/** To get Log2 of a FReal */
static float Log2(const float inValue){
return log2f(inValue);
}
static double Log2(const double inValue){
return log2(inValue);
}
/** To get Log of a FReal */
static float Log(const float inValue){
return logf(inValue);
}
static double Log(const double inValue){
return log(inValue);
}
/** To get atan2 of a 2 FReal, The return value is given in radians and is in the
range -pi to pi, inclusive. */
static float Atan2(const float inValue1,const float inValue2){
return atan2f(inValue1,inValue2);
}
static double Atan2(const double inValue1,const double inValue2){
return atan2(inValue1,inValue2);
}
/** To get Log2 of a FReal */
static float Log2(const float inValue){
return log2f(inValue);
}
static double Log2(const double inValue){
return log2(inValue);
}
/** To get sin of a FReal */
static float Sin(const float inValue){
return sinf(inValue);
}
static double Sin(const double inValue){
return sin(inValue);
}
/** To get atan2 of a 2 FReal, The return value is given in radians and is in the
range -pi to pi, inclusive. */
static float Atan2(const float inValue1,const float inValue2){
return atan2f(inValue1,inValue2);
}
static double Atan2(const double inValue1,const double inValue2){
return atan2(inValue1,inValue2);
}
/** To get asinf of a float. The result is in the range [0, pi]*/
static float ASin(const float inValue){
return asinf(inValue);
}
/** To get asinf of a double. The result is in the range [0, pi]*/
static double ASin(const double inValue){
return asin(inValue);
}
/** To get sin of a FReal */
static float Sin(const float inValue){
return sinf(inValue);
}
static double Sin(const double inValue){
return sin(inValue);
}
/** To get cos of a FReal */
static float Cos(const float inValue){
return cosf(inValue);
}
static double Cos(const double inValue){
return cos(inValue);
}
/** To get asinf of a float. The result is in the range [0, pi]*/
static float ASin(const float inValue){
return asinf(inValue);
}
/** To get asinf of a double. The result is in the range [0, pi]*/
static double ASin(const double inValue){
return asin(inValue);
}
/** To get arccos of a float. The result is in the range [0, pi]*/
static float ACos(const float inValue){
return acosf(inValue);
}
/** To get arccos of a double. The result is in the range [0, pi]*/
static double ACos(const double inValue){
return acos(inValue);
}
/** To get cos of a FReal */
static float Cos(const float inValue){
return cosf(inValue);
}
static double Cos(const double inValue){
return cos(inValue);
}
/** To get atan2 of a 2 FReal */
static float Fmod(const float inValue1,const float inValue2){
return fmodf(inValue1,inValue2);
}
/** return the floating-point remainder of inValue1 / inValue2 */
static double Fmod(const double inValue1,const double inValue2){
return fmod(inValue1,inValue2);
}
/** To get arccos of a float. The result is in the range [0, pi]*/
static float ACos(const float inValue){
return acosf(inValue);
}
/** To get arccos of a double. The result is in the range [0, pi]*/
static double ACos(const double inValue){
return acos(inValue);
}
/** To know if a variable is nan, based on the C++0x */
template <class TestClass>
static bool IsNan(const TestClass& value){
//volatile const TestClass* const pvalue = &value;
//return (*pvalue) != value;
return std::isnan(value);
}
/** To get atan2 of a 2 FReal */
static float Fmod(const float inValue1,const float inValue2){
return fmodf(inValue1,inValue2);
}
/** return the floating-point remainder of inValue1 / inValue2 */
static double Fmod(const double inValue1,const double inValue2){
return fmod(inValue1,inValue2);
}
/** To know if a variable is nan, based on the C++0x */
template <class TestClass>
static bool IsNan(const TestClass& value){
//volatile const TestClass* const pvalue = &value;
//return (*pvalue) != value;
return std::isnan(value);
}
/** To know if a variable is not inf, based on the C++0x */
template <class TestClass>
static bool IsFinite(const TestClass& value){
// return !(value <= std::numeric_limits<T>::min()) && !(std::numeric_limits<T>::max() <= value);
return std::isfinite(value);
}
/** To know if a variable is not inf, based on the C++0x */
template <class TestClass>
static bool IsFinite(const TestClass& value){
// return !(value <= std::numeric_limits<T>::min()) && !(std::numeric_limits<T>::max() <= value);
return std::isfinite(value);
}
/** A class to compute accuracy */
class FAccurater {
FReal l2Dot;
FReal l2Diff;
FReal max;
FReal maxDiff;
public:
FAccurater() : l2Dot(0), l2Diff(0), max(0), maxDiff(0) {
}
/** with inital values */
FAccurater(const FReal inGood[], const FReal inBad[], const int nbValues)
: l2Dot(0), l2Diff(0), max(0), maxDiff(0) {
add(inGood, inBad, nbValues);
}
/** Add value to the current list */
void add(const FReal inGood, const FReal inBad){
l2Diff += (inBad - inGood) * (inBad - inGood);
l2Dot += inGood * inGood;
/** A class to compute accuracy */
class FAccurater {
FReal l2Dot;
FReal l2Diff;
FReal max;
FReal maxDiff;
public:
FAccurater() : l2Dot(0), l2Diff(0), max(0), maxDiff(0) {
}
/** with inital values */
FAccurater(const FReal inGood[], const FReal inBad[], const int nbValues)
: l2Dot(0), l2Diff(0), max(0), maxDiff(0) {
add(inGood, inBad, nbValues);
}
/** Add value to the current list */
void add(const FReal inGood, const FReal inBad){
l2Diff += (inBad - inGood) * (inBad - inGood);
l2Dot += inGood * inGood;
max = Max(max , Abs(inGood));
maxDiff = Max(maxDiff, Abs(inGood-inBad));
}
/** Add array of values */
void add(const FReal inGood[], const FReal inBad[], const int nbValues){
for(int idx = 0 ; idx < nbValues ; ++idx){
add(inGood[idx],inBad[idx]);
}
}
/** Get the L2 norm */
FReal getL2Norm() const{
return Sqrt(l2Diff / l2Dot);
}
/** Get the inf norm */
FReal getInfNorm() const{
return maxDiff / max;
}
/** Print */
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FAccurater& inAccurater){
output << "[Error] L2Norm = " << inAccurater.getL2Norm() << " \t Inf = " << inAccurater.getInfNorm();
return output;
}
};
max = Max(max , Abs(inGood));
maxDiff = Max(maxDiff, Abs(inGood-inBad));
}
/** Add array of values */
void add(const FReal inGood[], const FReal inBad[], const int nbValues){
for(int idx = 0 ; idx < nbValues ; ++idx){
add(inGood[idx],inBad[idx]);
}
}
/** Get the L2 norm */
FReal getL2Norm() const{
return Sqrt(l2Diff / l2Dot);
}
/** Get the inf norm */
FReal getInfNorm() const{
return maxDiff / max;
}
/** Print */
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FAccurater& inAccurater){
output << "[Error] L2Norm = " << inAccurater.getL2Norm() << " \t Inf = " << inAccurater.getInfNorm();
return output;
}
void reset()
{
l2Dot = FReal(0);
l2Diff = FReal(0);;
max = FReal(0);
maxDiff = FReal(0);
}
};
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment