- 01
 - 02
 - 03
 - 04
 - 05
 - 06
 - 07
 - 08
 - 09
 - 10
 - 11
 - 12
 - 13
 - 14
 - 15
 - 16
 - 17
 - 18
 - 19
 - 20
 - 21
 - 22
 - 23
 - 24
 - 25
 - 26
 - 27
 - 28
 - 29
 - 30
 - 31
 - 32
 - 33
 - 34
 - 35
 - 36
 - 37
 - 38
 - 39
 - 40
 - 41
 - 42
 - 43
 - 44
 - 45
 - 46
 - 47
 - 48
 - 49
 - 50
 - 51
 - 52
 - 53
 - 54
 - 55
 - 56
 - 57
 - 58
 - 59
 - 60
 - 61
 - 62
 - 63
 - 64
 - 65
 - 66
 - 67
 - 68
 - 69
 - 70
 - 71
 - 72
 - 73
 - 74
 - 75
 - 76
 - 77
 - 78
 - 79
 - 80
 - 81
 
                        template < int Order, typename T >
struct PrefixSum {
  static inline void update ( T* a ) throw () { *a += *(a-1); PrefixSum < Order-1, T > :: update( a+1 ); }
};
template < typename T >
struct PrefixSum < 1,T > {
  static inline void update ( T* a ) throw () { *a += *(a-1); }
};
template < int P, int N, int Condition = 0 > 
struct Bpn {
  enum { value = N * ( Bpn < P-1, N-1, (N > P) > :: value - Bpn < P-1, N, (N > P-1) > :: value ) };
};
template < int P, int N > struct Bpn < P, N, !0 > { enum { value = 0 }; };
template < int P >        struct Bpn < P, 0,  0 > { enum { value = P ? 0 : 1 }; };
template < int P >        struct Bpn < P, 1,  0 > { enum { value = P & 1 ? 1 : -1 }; };
template < typename Ta, typename Tb, bool C > struct IfThenElse;
template < typename Ta, typename Tb > struct IfThenElse < Ta, Tb, true  > { typedef Ta TRes; };
template < typename Ta, typename Tb > struct IfThenElse < Ta, Tb, false > { typedef Tb TRes; };
template < int K, int I = 1 >
struct MomentSeries {
  typedef typename IfThenElse < MomentSeries<K,I+1>, MomentSeries<K,K>, (I<K) >::TRes SubT;
  static inline double accumulate ( double const* psum ) throw () {
    return Bpn < K, I > :: value * psum [ I + 1 ] + SubT :: accumulate ( psum );
  }
};
template < int K >
struct MomentSeries < K, K > {
  static inline double accumulate ( double const* psum ) throw () {
    return Bpn < K, K > :: value * psum [ K + 1 ];
  }
};
template < int Order >
struct MomentLoop {
  static inline void assign ( double *moments, size_t momentStride, double const* psum ) throw() {
    *(moments - momentStride) = MomentSeries < Order-1 > :: accumulate ( psum );
    MomentLoopAssign < Order-1 > :: assign ( moments - momentStride, momentStride, psum ); 
  }
};
template <>
struct MomentLoop < 1 > {
  static inline void assign ( double *moments, size_t momentStride, double const* psum ) throw() {
    moments [ 0 ] = MomentSeries < 1, 1 > :: accumulate ( psum );
    *(moments - momentStride) = psum [ 1 ];
  }
};
/**
 * Function computes a series of geometric moments by prefix summation method:
 * Zhou F., Kornerup P. Computing Moments by Prefix Sums. // VLSI Signal Proc. - 2000. - 25. - P. 5 - 17.
 * @param data is first data elemet address.
 * @param ndataItems is number of data items.
 * @param dataStride is the number of elements between two neighbor data items, 
 *        in case of simple array dataStride is equal to 1.
 * @param moments is address of 0-order moment.
 * @param momentStride is number of elements between two neigbor moment items,
 *        in case of consequtive moments placement, momentStride is equal to 1.
 */
template < int Order, class OutPolicy >
inline void psmoment ( double const* data, 
                       size_t const  ndataItems,
                       size_t const  dataStride,
                       double*       moments,
                       size_t const  momentStride ) throw() { 
  double psum [ Order+1 ] = { 0 };
  // Initialize prefix sum
  for ( size_t i = ndataItems, j = ndataItems * dataStride; i; ) {
    --i; j -= dataStride; psum [ 0 ] = data [ j ];
    PrefixSum < Order, double > :: update ( psum+1 );
  }
  // convert psum to moment values 
  OutPolicy :: assign ( moments + Order * momentStride, momentStride, psum );
}