/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/

/*! \file simd_type.h

    \brief SIMD type derived from std::simd

    To use this header, an implementation of std::simd has to be
    installed, and the -std=c++17 option is needed as well.
    It has been tried with clang++ and g++; you'll need a
    recent version.

*/

#ifndef VSPLINE_STD_SIMD_TYPE_H
#define VSPLINE_STD_SIMD_TYPE_H

#ifndef VSPLINE_VECTOR_NBYTES

#define VSPLINE_VECTOR_NBYTES 64

#endif

#include <iostream>
#include <experimental/simd>

namespace vspline
{

/// class template simd_type provides a fixed-size SIMD type.
/// This implementation of vspline::simd_type uses std::simd as
/// base class. This class is used as a stand-in for Vc::SimdArray
/// - it does not cover the whole interface, but a reasonably
/// large subset - my choice of SimdArray functionality is what
/// I need in lux/vspline.
/// Most of the 'loop variant' of simd_type has been ported to
/// use std::simd instead, making use of std::simd's
/// - constructors
/// - copy_from and copy_to
/// - masks and where expressions
/// - operator functions
/// - overloads for several mathematical functions
/// - min/max

template < typename _value_type ,
           std::size_t _vsize >
struct simd_type
: private std::experimental::simd
          < _value_type ,
            std::experimental::simd_abi::fixed_size < _vsize >
          >
{
  typedef std::experimental::simd_abi::fixed_size < _vsize > abi_t ;

  typedef std::size_t size_type ;
  typedef _value_type value_type ;
  static const size_type vsize = _vsize ;

  typedef std::experimental::simd < value_type , abi_t > base_t ;
  typedef std::experimental::simd < int , abi_t > index_type ;
  using typename base_t::mask_type ;

  // provide the size as a constexpr

  static constexpr size_type size()
  {
    return vsize ;
  }

  typedef mask_type MaskType ;
  typedef index_type IndexType ;

  // operator[] is mapped to std::simd element access

  using base_t::operator[] ;

  // assignment from a value_type. The assignment is coded as a loop,
  // but it should be obvious to the compiler's loop vectorizer that
  // the loop is a 'SIMD operation in disguise', so here we have the
  // first appearance of 'goading'.

  simd_type & operator= ( const value_type & rhs )
  {
    to_base() = base_t ( rhs ) ;
    return *this ;
  }

  // c'tor from value_type. We use the assignment operator for
  // initialization.

  simd_type ( const value_type & ini )
  : base_t ( ini )
  { }

  simd_type ( const base_t & ini )
  : base_t ( ini )
  { }

  // these two c'tors are left in default mode

  simd_type() = default ;
  simd_type ( const simd_type & ) = default ;

  base_t & to_base()
  {
    return ( * static_cast < base_t * > ( this ) ) ;
  }

  const base_t & to_base() const
  {
    return ( * static_cast < const base_t * const > ( this ) ) ;
  }

  // assignment from equally-sized container. Most containers use std::size_t
  // for the template argument defining the number of elements they hold,
  // but some (notably vigra::TinyVector) use int, which is probably a relic
  // from times when non-type template arguments were of a restricted type
  // set only. By providing a specialization for SIZE_TYPE int, we make
  // equally-sized vigra::TinyVectors permitted initializers.
  // the c'tor from an equally-sized container also uses the corresponding
  // operator= overload, so we use one macro for both.
  // we also need two different variants of vsize for g++; clang++ accepts
  // size_type vsize for both places where VSZ is used, but g++ requires
  // an integer.
  // Note that the rhs can use any elementary type which can be legally
  // assigned to value_type. This allows transport of information from
  // differently typed objects, but there are no further constraints on
  // the types involved, which may degrade precision. It's the user's
  // responsibility to make sure such assignments have the desired effect
  // and overload them if necessary.

  #define BUILD_FROM_CONTAINER(SIZE_TYPE,VSZ) \
    template < typename U , template < typename , SIZE_TYPE > class V > \
    simd_type & operator= ( const V < U , VSZ > & rhs ) \
    { \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        to_base() [ i ] = value_type ( rhs [ i ] ) ; \
      return *this ; \
    } \
    template < typename U , template < typename , SIZE_TYPE > class V > \
    simd_type ( const V < U , VSZ > & ini ) \
    { \
      *this = ini ; \
    }

  BUILD_FROM_CONTAINER(std::size_t,vsize)

  #undef BUILD_FROM_CONTAINER

  static const simd_type iota()
  {
    simd_type result ;
    for ( size_type i = 0 ; i < vsize ; i++ )
      result [ i ] = value_type ( i ) ;
    return result ;
  }

  // mimick Vc's IndexesFromZero. This function produces an index
  // vector filled with indexes starting with zero.

  static const index_type IndexesFromZero()
  {
    typedef typename index_type::value_type IT ;
    static const IT ceiling = std::numeric_limits < IT > :: max() ;
    assert ( ( vsize - 1 ) <= std::size_t ( ceiling ) ) ;

    index_type ix ;
    for ( size_type i = 0 ; i < vsize ; i++ )
      ix [ i ] = int ( i ) ;
    return ix ;
  }

  // variant which starts from a different starting point and optionally
  // uses steps other than one.

  static const index_type IndexesFrom ( std::size_t start ,
                                        std::size_t step = 1 )
  {
    typedef typename index_type::value_type IT ;
    static const IT ceiling = std::numeric_limits < IT > :: max() ;
    assert ( start + ( vsize - 1 ) * step <= std::size_t ( ceiling ) ) ;

    return ( IndexesFromZero() * int(step) ) + int(start) ;
  }

  // functions Zero and One produce simd_type objects filled with
  // 0, or 1, respectively

  static const simd_type Zero()
  {
    return simd_type ( value_type ( 0 ) ) ;
  }

  static const simd_type One()
  {
    return simd_type ( value_type ( 1 ) ) ;
  }

  // echo the vector to a std::ostream, read it from an istream

  friend std::ostream & operator<< ( std::ostream & osr ,
                                     simd_type it )
  {
    osr << "(" ;
    for ( size_type i = 0 ; i < vsize - 1 ; i++ )
      osr << it [ i ] << ", " ;
    osr << it [ vsize - 1 ] << ")" ;
    return osr ;
  }

  friend std::istream & operator>> ( std::istream & isr ,
                                     simd_type it )
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      isr >> it [ i ] ;
    return isr ;
  }

  // memory access functions, which load and store vector data.
  // We start out with functions transporting data from memory into
  // the simd_type. Some of these operations have corresponding
  // c'tors which use the member function to initialize to_base().

  // load delegates to std::simd::copy_from. TODO: consider
  // overalignment

  void load ( const value_type * const p_src )
  {
    to_base().copy_from ( p_src ,
                          std::experimental::element_aligned_tag() ) ;
  }

  // std::simd does not offer gather/scatter, but it offers a
  // c'tor taking a functor to set the elements. In theory, this
  // is a good idea, because the optimizer might realize that
  // the sum of invocations of the functor can be represented
  // by a gather/scatter operation, but how well this works is
  // a different matter and mileage varies.

#define GS_LAMBDA

#ifdef GS_LAMBDA

  template < typename index_type >
  void gather ( const value_type * const p_src ,
                const index_type & indexes )
  {
    // assign base_t object created by gen-type c'tor
    to_base() = base_t ( [&] ( const size_t & i )
                    { return p_src [ indexes[i] ] ; } ) ;
  }

#else

  template < typename index_type >
  void gather ( const value_type * const p_src ,
                const index_type & indexes )
  {
    for ( std::size_t i = 0 ; i < vsize ; i++ )
      (*this)[i] = p_src [ indexes [ i ] ] ;
  }

#endif

  // c'tor from pointer and indexes, uses gather

  template < typename index_type >
  simd_type ( const value_type * const p_src ,
              const index_type & indexes )
  {
    gather ( p_src , indexes ) ;
  }

  // store saves the content of the vector to memory

  void store ( value_type * const p_trg ) const
  {
    to_base().copy_to ( p_trg ,
                        std::experimental::element_aligned_tag() ) ;
  }

  // scatter is the reverse operation to gather

#ifdef GS_LAMBDA

  template < typename index_type >
  void scatter ( value_type * const p_trg ,
                 const index_type & indexes ) const
  {
    // gen-type c'tor is only used for side effects; let the compiler
    // figure out that the result is unused.

    base_t dummy ( [&] ( const size_t & i )
                   { return p_trg [ indexes[i] ] = to_base()[i] ; } ) ;
  }

#else

  template < typename index_type >
  void scatter ( value_type * const p_trg ,
                 const index_type & indexes ) const
  {
    for ( std::size_t i = 0 ; i < vsize ; i++ )
      p_trg [ indexes [ i ] ] = (*this)[i] ;
  }

#endif

  // 'regular' gather and scatter, accessing strided memory so that the
  // first address visited is p_src/p_trg, and successive addresses are
  // 'step' apart - in units of T. Might also be done with goading, the
  // loop should autovectorize.

  void rgather ( const value_type * const p_src ,
                 const std::size_t & step )
  {
    auto indexes = IndexesFrom ( 0 , step ) ;
    gather ( p_src , indexes ) ;
  }

  void rscatter ( value_type * p_trg ,
                  const std::size_t & step ) const
  {
    auto indexes = IndexesFrom ( 0 , step ) ;
    scatter ( p_trg , indexes ) ;
  }

  // apply functions from namespace std to each element in a vector,
  // or to each corresponding set of elements in a set of vectors
  // - going up to three for fma.
  // many standard functions autovectorize well. Note that the
  // autovectorization of standard functions often needs additional
  // compiler flags, like, e.g., -fno-math-errno for clang++, to
  // produce hardware SIMD instructions.

  #define BROADCAST_STD_FUNC(FUNC) \
    friend simd_type FUNC ( simd_type arg ) \
    { \
      return FUNC ( arg.to_base() ) ; \
    }

    // TODO: getting zero back for negative args, hence no BROADCAST_STD_FUNC
    // this happens with clang++ only, I opened an issue with VcDevel/std-simd:
    // https://github.com/VcDevel/std-simd/issues/31

//   BROADCAST_STD_FUNC(abs)

  friend simd_type abs ( simd_type arg )
  {
    arg ( arg < 0 ) = - arg ;
    return arg ;
  }

  BROADCAST_STD_FUNC(trunc)

  BROADCAST_STD_FUNC(round)
  BROADCAST_STD_FUNC(floor)
  BROADCAST_STD_FUNC(ceil)
  BROADCAST_STD_FUNC(log)
  BROADCAST_STD_FUNC(exp)
  BROADCAST_STD_FUNC(sqrt)

  // the support for autovectorization of trigonometric functions is
  // sketchy - e.g. the clang++ reference does not mention them as
  // functions which autovectorize. Vc offers hand-coded trigonometric
  // functions which might be worth while porting to std::simd, but so
  // far I haven't seen this happen. In my application, this results in
  // bad performance when these functions are used.

  BROADCAST_STD_FUNC(tan)
  BROADCAST_STD_FUNC(asin)
  BROADCAST_STD_FUNC(acos)
  BROADCAST_STD_FUNC(atan)

  // TODO: odd: with clang++, sin and cos don't perform as expected;
  // using a loop does the trick:

#ifdef __clang__

  friend simd_type cos ( simd_type arg )
  {
    simd_type result ;
    for ( std::size_t i = 0 ; i < size() ; i++ )
      result[i] = std::cos ( arg[i] ) ;
    return result ;
  }

  friend simd_type sin ( simd_type arg )
  {
    simd_type result ;
    for ( std::size_t i = 0 ; i < size() ; i++ )
      result[i] = std::sin ( arg[i] ) ;
    return result ;
  }

#else

  BROADCAST_STD_FUNC(sin)
  BROADCAST_STD_FUNC(cos)

#endif

  #undef BROADCAST_STD_FUNC

  #define BROADCAST_STD_FUNC2(FUNC) \
    friend simd_type FUNC ( simd_type arg1 , \
                            simd_type arg2 ) \
    { \
      return FUNC ( arg1.to_base() , arg2.to_base() ) ; \
    }

  // a short note on atan2: Vc provides a hand-written vectorized version
  // of atan2 which is especially fast and superior to autovectorized code.

  BROADCAST_STD_FUNC2(atan2)
  BROADCAST_STD_FUNC2(pow)

  #undef BROADCAST_STD_FUNC2

  #define BROADCAST_STD_FUNC3(FUNC) \
    friend simd_type FUNC ( simd_type arg1 , \
                            simd_type arg2 , \
                            simd_type arg3 ) \
    { \
      return FUNC ( arg1.to_base() , arg2.to_base() , arg3.to_base() ) ; \
    }

  BROADCAST_STD_FUNC3(fma)

  #undef BROADCAST_STD_FUNC3

  // macros used for the parameter 'CONSTRAINT' in the definitions
  // further down. Some operations are only allowed for integral types
  // or boolans. This might be enforced by enable_if, here we use a
  // static_assert with a clear error message.
  // TODO: might relax constraints by using 'std::is_convertible'

  #define INTEGRAL_ONLY \
    static_assert ( std::is_integral < value_type > :: value , \
                    "this operation is only allowed for integral types" ) ;

  #define BOOL_ONLY \
    static_assert ( std::is_same < value_type , bool > :: value , \
                    "this operation is only allowed for booleans" ) ;

  // augmented assignment operators. Some operators are only applicable
  // to specific data types, which is enforced by 'CONSTRAINT'.
  // One might consider widening the scope by making these operator
  // functions templates and accepting arbitrary indexable types.
  // Only value_type and simd_type itto_base() are taken as rhs arguments.

  #define OPEQ_FUNC(OPFUNC,OPEQ,CONSTRAINT) \
    simd_type & OPFUNC ( value_type rhs ) \
    { \
      CONSTRAINT \
      to_base() OPEQ rhs ; \
      return *this ; \
    } \
    simd_type & OPFUNC ( simd_type rhs ) \
    { \
      CONSTRAINT \
      to_base() OPEQ rhs.to_base() ; \
      return *this ; \
    }

  OPEQ_FUNC(operator+=,+=,)
  OPEQ_FUNC(operator-=,-=,)
  OPEQ_FUNC(operator*=,*=,)
  OPEQ_FUNC(operator/=,/=,)

  OPEQ_FUNC(operator%=,%=,INTEGRAL_ONLY)
  OPEQ_FUNC(operator&=,&=,INTEGRAL_ONLY)
  OPEQ_FUNC(operator|=,|=,INTEGRAL_ONLY)
  OPEQ_FUNC(operator^=,^=,INTEGRAL_ONLY)
  OPEQ_FUNC(operator<<=,<<=,INTEGRAL_ONLY)
  OPEQ_FUNC(operator>>=,>>=,INTEGRAL_ONLY)

  #undef OPEQ_FUNC

  // we use a simple scheme for type promotion: the promoted type
  // of two values should be the same as the type we would receive
  // when adding the two values. That's standard C semantics, but
  // it won't widen the result type to avoid overflow or increase
  // precision - such conversions have to be made by user code if
  // necessary.

  #define C_PROMOTE(A,B)  \
  typename std::conditional \
            < std::is_same < A , B > :: value , \
              A , \
              decltype (   std::declval < A > () \
                          + std::declval < B > () ) \
            > :: type

  // binary operators and left and right scalar operations with
  // value_type, unary operators -, ! and ~

  #define OP_FUNC(OPFUNC,OP,CONSTRAINT) \
  template < typename RHST , \
             typename = typename std::enable_if \
                       < std::is_fundamental < RHST > :: value \
                       > :: type \
           > \
  simd_type < C_PROMOTE ( value_type , RHST ) , size() > \
  OPFUNC ( simd_type < RHST , vsize > _rhs ) const \
  { \
    CONSTRAINT \
    simd_type < C_PROMOTE ( value_type , RHST ) , vsize > lhs ( *this ) ; \
    simd_type < C_PROMOTE ( value_type , RHST ) , vsize > rhs ( _rhs ) ; \
    return lhs.to_base() OP rhs.to_base() ; \
  } \
  template < typename RHST , \
             typename = typename std::enable_if \
                       < std::is_fundamental < RHST > :: value \
                       > :: type \
           > \
  simd_type < C_PROMOTE ( value_type , RHST ) , vsize > \
  OPFUNC ( RHST _rhs ) const \
  { \
    CONSTRAINT \
    simd_type < C_PROMOTE ( value_type , RHST ) , vsize > lhs ( *this ) ; \
    C_PROMOTE ( value_type , RHST ) rhs ( _rhs ) ; \
    return lhs.to_base() OP rhs ; \
  } \
  template < typename LHST , \
             typename = typename std::enable_if \
                       < std::is_fundamental < LHST > :: value \
                       > :: type \
           > \
  friend simd_type < C_PROMOTE ( LHST , value_type ) , vsize > \
  OPFUNC ( LHST _lhs , simd_type _rhs ) \
  { \
    CONSTRAINT \
    C_PROMOTE ( value_type , LHST ) lhs ( _lhs ) ; \
    simd_type < C_PROMOTE ( LHST , value_type ) , vsize > rhs ( _rhs ) ; \
    return lhs OP rhs.to_base() ; \
  }

  // binary operators and left and right scalar operations with
  // value_type, unary operators -, ! and ~

  // #define OP_FUNC(OPFUNC,OP,CONSTRAINT) \
  //   simd_type OPFUNC ( simd_type rhs ) const \
  //   { \
  //     CONSTRAINT \
  //     return this->to_base() OP rhs.to_base() ; \
  //   } \
  //   simd_type OPFUNC ( value_type rhs ) const \
  //   { \
  //     CONSTRAINT \
  //     return this->to_base() OP rhs ; \
  //   } \
  //   friend simd_type OPFUNC ( value_type lhs , \
  //                             simd_type rhs ) \
  //   { \
  //     CONSTRAINT \
  //     return lhs OP rhs.to_base() ; \
  //   }
  //
  OP_FUNC(operator+,+,)
  OP_FUNC(operator-,-,)
  OP_FUNC(operator*,*,)
  OP_FUNC(operator/,/,)

  OP_FUNC(operator%,%,INTEGRAL_ONLY)
  OP_FUNC(operator&,&,INTEGRAL_ONLY)
  OP_FUNC(operator|,|,INTEGRAL_ONLY)
  OP_FUNC(operator^,^,INTEGRAL_ONLY)
  OP_FUNC(operator<<,<<,INTEGRAL_ONLY)
  OP_FUNC(operator>>,>>,INTEGRAL_ONLY)

  OP_FUNC(operator&&,&&,BOOL_ONLY)
  OP_FUNC(operator||,||,BOOL_ONLY)

  #undef OP_FUNC

  #define OP_FUNC(OPFUNC,OP,CONSTRAINT) \
    simd_type OPFUNC() const \
    { \
      return OP this->to_base() ; \
    }

  OP_FUNC(operator-,-,)
  OP_FUNC(operator!,!,BOOL_ONLY)
  OP_FUNC(operator~,~,INTEGRAL_ONLY)

  #undef OP_FUNC

  // provide methods to produce a mask on comparing a vector
  // with another vector or a value_type.

  #define COMPARE_FUNC(OP,OPFUNC) \
  friend mask_type OPFUNC ( simd_type lhs , \
                            simd_type rhs ) \
  { \
    return lhs.to_base() OP rhs.to_base() ; \
  } \
  friend mask_type OPFUNC ( simd_type lhs , \
                            value_type rhs ) \
  { \
    return lhs.to_base() OP rhs ; \
  } \
  friend mask_type OPFUNC ( value_type lhs , \
                            simd_type rhs ) \
  { \
    return lhs OP rhs.to_base() ; \
  }

  COMPARE_FUNC(<,operator<) ;
  COMPARE_FUNC(<=,operator<=) ;
  COMPARE_FUNC(>,operator>) ;
  COMPARE_FUNC(>=,operator>=) ;
  COMPARE_FUNC(==,operator==) ;
  COMPARE_FUNC(!=,operator!=) ;

  #undef COMPARE_FUNC

  // note: std::simd's mask_type has associated functions any_of, all_of
  // and none_of, so we needn't define them for this backend

  // next we define a masked vector as an object holding two references:
  // one reference to a mask type, determining which of the vector's
  // elements will be 'open' to an effect, and one reference to a vector,
  // which will be affected by the operation.
  // The resulting object will only be viable as long as the referred-to
  // mask and vector are alive - it's meant as a construct to be processed
  // in the same scope, as the lhs of an assignment, typically using
  // notation introduced by Vc: a vector's operator() is overloaded to
  // to produce a masked_type when called with a mask_type object, and
  // the resulting masked_type object is then assigned to.
  // Note that this does not have any effect on those values in 'whither'
  // for which the mask is false. They remain unchanged.

  typedef std::experimental::where_expression < mask_type , base_t > we_t ;

  struct masked_type
  {
    mask_type whether ;   // if the mask is true at whether[i]
    simd_type & whither ; // whither[i] will be assigned to

    masked_type ( mask_type _whether ,
                  simd_type & _whither )
    : whether ( _whether ) ,
      whither ( _whither )
      { }

    // for the masked vector, we define the complete set of assignments:

    #define OPEQ_FUNC(OPFUNC,OPEQ,CONSTRAINT) \
      simd_type & OPFUNC ( value_type rhs ) \
      { \
        CONSTRAINT \
        we_t ( whether , whither.to_base() ) OPEQ rhs ; \
        return whither ; \
      } \
      simd_type & OPFUNC ( simd_type rhs ) \
      { \
        CONSTRAINT \
        we_t ( whether , whither.to_base() ) OPEQ rhs.to_base() ; \
        return whither ; \
      }

    OPEQ_FUNC(operator=,=,)
    OPEQ_FUNC(operator+=,+=,)
    OPEQ_FUNC(operator-=,-=,)
    OPEQ_FUNC(operator*=,*=,)
    OPEQ_FUNC(operator/=,/=,)
    OPEQ_FUNC(operator%=,%=,INTEGRAL_ONLY)
    OPEQ_FUNC(operator&=,&=,INTEGRAL_ONLY)
    OPEQ_FUNC(operator|=,|=,INTEGRAL_ONLY)
    OPEQ_FUNC(operator^=,^=,INTEGRAL_ONLY)
    OPEQ_FUNC(operator<<=,<<=,INTEGRAL_ONLY)
    OPEQ_FUNC(operator>>=,>>=,INTEGRAL_ONLY)

  #undef OPEQ_FUNC

  #undef INTEGRAL_ONLY
  #undef BOOL_ONLY
  
  } ;

  // mimicking Vc, we define operator() with a mask_type argument
  // to produce a masked_type object, which can be used later on to
  // masked-assign to the referred-to vector. With this definition
  // we can use the same syntax Vc uses, e.g. v1 ( v1 > v2 ) = v3
  // This helps write code which compiles with Vc and without,
  // because this idiom is 'very Vc'.

  masked_type operator() ( mask_type mask )
  {
    return masked_type ( mask , *this ) ;
  }

  // member functions at_least and at_most. These functions provide the
  // same functionality as max, or min, respectively. Given simd_type X
  // and some threshold Y, X.at_least ( Y ) == max ( X , Y )
  // Having the functionality as a member function makes it easy to
  // implement, e.g., min as: min ( X , Y ) { return X.at_most ( Y ) ; }

  #define CLAMP(FNAME,REL) \
    simd_type FNAME ( simd_type threshold ) const \
    { \
      return REL ( to_base() , threshold.to_base() ) ; \
    } \
    simd_type FNAME ( value_type threshold ) const \
    { \
      return REL ( to_base() , threshold ) ; \
    } \

  CLAMP(at_least,max)
  CLAMP(at_most,min)

  #undef CLAMP

  // sum of vector elements. Note that there is no type promotion; the
  // summation is done to value_type. Caller must make sure that overflow
  // is not a problem.

  value_type sum() const
  {
    value_type s ( 0 ) ;
    for ( std::size_t e = 0 ; e < vsize ; e++ )
      s += (*this) [ e ] ;
    return s ;
  }
} ;

} ;

#endif // #define VSPLINE_SIMD_TYPE_H