toolbox/arrayOp.hh

Go to the documentation of this file.

#ifndef ArrayOp_hh
#define ArrayOp_hh

#include <stdarg.h>
#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;
  }

  // ******************************************************************* 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;
  }

} // namespace std

namespace concepts {

  // ************************************************************* makeArray **

  template<class F>
  Array<F> makeArray(uint n, const F& first, ...) {
    Array<F> data(n);
    data[0] = first;

    va_list param;
    va_start(param, first);
    for(uint i = 1; i < n; ++i)
      data[i] = va_arg(param, F);
    va_end(param);  
    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));
  }

} // namespace concepts

#endif // ArrayOp_hh