n

Src/Kernels/Rotation/FRotationOriginalKernel.hpp

@@ -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");

Src/Utils/FComplexe.hpp

@@ -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.

Src/Utils/FMath.hpp

@@ -23,185 +23,203 @@
 #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 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 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);
+  }
 
 
+  /** 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;
@@ -253,6 +271,15 @@ struct FMath{
             return output;
         }
     };
+
+    void reset()
+    {
+      l2Dot = FReal(0);
+      l2Diff = FReal(0);;
+      max = FReal(0);
+      maxDiff = FReal(0);
+    }
+  };
 };