- 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 );
}