arrayOp.hh

Go to the documentation of this file.

 
#ifndef ArrayOp_hh
#define ArrayOp_hh
 
#include <stdarg.h>
#include <initializer_list>
#include "toolbox/array.hh"
 
// debugging
#define ArrayProduct_D 0
 
namespace std {
 
  // *************************************************************** product **
 
  template<typename F>
  F product(const concepts::Array<F>& a) {
    const F* d = (const F*)a;
    F val = F(1);
    for(uint i = a.size(); i--; ) val *= *d++;
    return val;
  }
 
  template<typename F>
  F product(const concepts::Array<F>& a, uint j) {
    const F* d = (const F*)a;
    F val = F(1);
    for(uint i = 0; i < a.size(); ++i, ++d)
      if (i != j) val *= *d;
    return val;
  }
 
  template<typename F>
  F product(const concepts::Array<F>& a, uint j, uint k) {
    const F* d = (const F*)a;
    F val = F(1);
    for(uint i = 0; i < a.size(); ++i, ++d) 
      if (i != j && i != k) val *= *d;
    return val;
  }
 
  // ******************************************************************* pow **
 
  template<typename F>
  concepts::Array<F> pow(const concepts::Array<F>& a, const F e) {
    concepts::Array<F> res(a.size());
    const F* d = a; F* o = res;
    for(uint i = a.size(); i--; ) *o++ = std::pow(*d++, e);
    return res;
  }
 
  // ******************************************************************* abs **
 
  template<typename F>
  concepts::Array<F> abs(const concepts::Array<F>& a) {
    concepts::Array<F> res(a.size());
    const F* d = a; F* e = res;
    for(uint i = a.size(); i--; ) *e++ = std::abs(*d++);
    return res;
  }
 
  // ******************************************************************* min **
 
  template<typename F>
  F min(const concepts::Array<F>& a) {
    conceptsAssert(a.size() > 0, concepts::Assertion());
    const F* d = a;
    F val = *d++;
    for(uint i = a.size()-1; i--; ) 
      val = std::min(val, *d++);
    return val;
  }
 
  // ******************************************************************* max **
 
  template<typename F>
  F max(const concepts::Array<F>& a) {
    conceptsAssert(a.size() > 0, concepts::Assertion());
    const F* d = a;
    F val = *d++;
    for(uint i = a.size()-1; i--; )
      val = std::max(val, *d++);
    return val;
  }
 
} // namespace std
 
namespace concepts {
 
  // ************************************************************* makeArray **
 
  template <class F>
  Array<F> makeArray(std::initializer_list<F> list)
  {
    uint n = list.size();
    Array<F> data(n);
    uint idx = 0;
    for( F elem : list )
      data[idx++] = elem;
    return data;
  }
 
  // ******************************************************* chebychevPoints **
 
  inline void chebychevPoints(concepts::Array<Real>& p) {
    uint n = p.size();
    const Real f = M_PI / (2*n);
    Real *x = (Real*)p;
    while(n) *x++ = cos(f*(2*--n+1));
  }
 
  // ******************************************************** componentArray **
 
  template <class F, uint dim>
  Array<F> componentArray(const Array<Point<F,dim>>& a, uint i) {
    conceptsAssert(i < dim, concepts::Assertion());
    uint n = a.size();
    Array<F> data(n);
    for(uint j = 0; j < n; ++j)
      data[j] = a[j][i];
    return data;
  }
 
} // namespace concepts
 
#endif // ArrayOp_hh