simd_type.h

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    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 using small loops
 
    vspline can use Vc for explicit vectorization, and at the time of
    this writing, this is usually the best option. But Vc is not available
    everywhere, or it's use may be unwanted. To help with such situations,
    vspline defines it's own 'SIMD' type, which is implemented as a simple
    C vector and small loops operating on it. If these constructs are
    compiled with compilers capable of autovectorization (and with the
    relevent flags activating use of SIMD instruction sets like AVX)
    the resulting code will oftentimes be 'proper' SIMD code, because
    the small loops are presented so that the compiler can easily recognize
    them as potential clients of loop vectorization. I call this technique
    'goading': By presenting the data flow in deliberately vector-friendly
    format, the compiler is more likely to 'get it'.
 
    class template simd_type is designed to provide an interface similar
    to Vc::SimdArray, to be able to use it as a drop-in replacement.
    It aims to provide those SIMD capabilities which are actually used by
    vspline and is *not* a complete replacement for Vc::SimdArray.
 
    Wherever possible, the code is as simple as possible, avoiding frills
    and trickery which might keep the compiler from recognizing potentially
    auto-vectorizable constructs. The resulting code is - in my limited
    experience - often not too far from explicit SIMD code. Some constructs
    do actually produce binary which is en par with code using Vc, namely
    such code which does not use gather, scatter or masked operations.
    So b-spline prefiltering, restoration of original data, and general
    filtering is very fast, while code involving b-spline evaluation
    shows a speed penalty, since vectorized b-spline evaluation (as coded
    in vspline) relies massively on gather operations of a kind which seem
    not to be auto-vectorized into binary gather commands - this is my guess,
    I have not investigated the binary closely.
 
    The code presented here adds some memory access functions which are
    not present in Vc::SimdArray, namely strided load/store operations
    and load/store using functors.
 
    Note that I use clang++ most of the time, and the code has evolved to
    produce fast binary with clang++. Your mileage will vary with other
    compilers.
 
    Class vspline::simd_type is actually quite similar to vigra::TinyVector
    which also stores in a plain C array and provides arithmetic. But that
    type is quite complex, using CRTP with a base class, explicitly coding
    loop unrolling, catering for deficient compilers and using vigra's
    sophisticated type promotion mechanism. vspline::simd_type on the other
    hand is stripped down to the bare essentials, to make the code as simple
    as possible, in the hope that 'goading' will indeed work. It replaces
    vspline's previous SIMD type, vspline::simd_tv, which was derived
    from vigra::TinyVector.
 
    One word of warning: the lack of type promotion requires you to pick
    a value_type of sufficient precision and capacity for the intended
    operation. In other words: you won't get an int when multiplying two
    shorts.
 
    Note also that this type is intended for *horizontal* vectorization,
    and you'll get the best results when picking a vector size which is
    a small-ish power of two - preferably at least the number of values
    of the given value_type which a register of the intended vector ISA
    will contain.
 
    vspline uses TinyVectors of SIMD data types, but their operations are
    coded with loops over the TinyVector's elements throughout vspline's
    code base. In vspline's opt directory, you can find 'xel_of_vector.h',
    which can provide overloads for all operator functions involving
    TinyVectors of vspline::simd_type - or, more generally, small
    aggregates of vector data. Please see this header's comments for
    more detailed information.
 
    Note also that throughout vspline, there is almost no explicit use of
    vspline::simd_type. vspline picks appropriate SIMD data types with
    mechanisms 'one level up', coded in vector.h. vector.h checks if use
    of Vc is possible and whether Vc can vectorize a given type, and
    produces a 'simdized type', which you mustn't confuse with a simd_type.
*/
 
#ifndef VSPLINE_SIMD_TYPE_H
#define VSPLINE_SIMD_TYPE_H
 
#include <iostream>
 
namespace vspline
{
 
template < typename T >
using is_scalar = typename std::integral_constant < bool ,
     std::is_fundamental < T > ::value
  || std::is_same < T , bool > :: value > :: type ;
 
/// class template simd_type provides a fixed-size container type for
/// small sets of fundamentals which are stored in a C vector.
/// The type offers arithmetic capabilities which are implemented by
/// using loops over the elements in the vector, expecting that the
/// compiler will autovectorize these loops into 'proper' SIMD code.
/// The interface of this type is modelled to be compatible with
/// Vc's SimdArray. Unfortunately, Vc::SimdArray requires additional
/// template arguments, so at times it's difficult to use the two types
/// instead of each other. The interface compatibility does *not* mean
/// that the arithmetic will produce the same results - this is intended
/// but neither tested nor enforced.
 
template < typename _value_type ,
           std::size_t _vsize >
class simd_type
{
  // storage of data is in a simple C array. This array is private,
  // and the only access to it is via member functions. Using a plain
  // C array will not overalign the data, even though this may be
  // desirable, especially on older hardware. We rely on the
  // compiler to handle this as efficiently as possible.
 
  _value_type _store [ _vsize ] ;
 
public:
 
  // some common typedefs for a container type
 
  typedef std::size_t size_type ;
  typedef _value_type value_type ;
  static const size_type vsize = _vsize ;
  static const int ivsize = _vsize ;      // finessing for g++
 
  // provide the size as a constexpr
 
  static constexpr size_type size()
  {
    return vsize ;
  }
 
  // types used for masks and index vectors. In terms of 'true' SIMD
  // arithmetics, these definitions may not be optimal - especially the
  // definition of a mask as a simd_type of bool is questionable - one
  // might consider using a bit field or a sufficiently large integral
  // type. But using a simd_type of bool makes processing simple, in
  // a way it's the 'generic' mask type, whereas SIMD masks used by
  // the hardware are the truly 'exotic' types. The problem here is
  // the way C++ encodes booleans - they are usually encoded as some
  // smallish integral type, rather than a single bit.
 
  // we define both the 'old school' and the 'camel case' variants
 
  typedef simd_type < bool , vsize > mask_type ;
  typedef simd_type < int , vsize > index_type ;
 
  typedef simd_type < bool , vsize > MaskType ;
  typedef simd_type < int , vsize > IndexType ;
 
  // operator[] is mapped to ordinary element access on a C vector.
  // This is the only place where we explicitly access _store, the
  // remainder of the code uses operator[].
 
  value_type & operator[] ( const size_type & i )
  {
    return _store[i] ;
  }
 
  const value_type & operator[] ( const size_type & i ) const
  {
    return _store[i] ;
  }
 
  // 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 )
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      (*this) [ i ] = rhs ;
    return *this ;
  }
 
  // c'tor from value_type. We use the assignment operator for
  // initialization.
 
  simd_type ( const value_type & ini )
  {
    *this = ini ;
  }
 
  // these two c'tors are left in default mode
 
  simd_type() = default ;
  simd_type ( const simd_type & ) = default ;
 
  // 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++ ) \
        (*this) [ i ] = 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)
  BUILD_FROM_CONTAINER(int,ivsize)
 
  #undef BUILD_FROM_CONTAINER
 
  // this operator= overload caters for Vc::SimdArray as rhs type.
  // Vc::SimdArray requires additional template arguments - 'an
  // unfortunate implementation detail shining through'.
  // We don't want to refer explicitly to Vc::SimdArray here, so
  // we just code a template which will 'catch' it. In result, we
  // can now easily assign a Vc::SimdArray to a vspline::simd_type.
  // The reverse operation can't be coded like this, though, since
  // assignment operators and c'tors have to be coded as members of
  // the 'receiving end'. See 'offload' for an efficient way of
  // moving simd_type data to Vc::SimdArrays
 
  template < typename U ,
             typename V ,
             std::size_t Wt ,
             template < typename , std::size_t ,
                        typename , std::size_t > class W >
  simd_type & operator= ( const W < U , vsize , V , Wt > & rhs )
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      (*this) [ i ] = rhs [ i ] ; 
    return *this ;
  }
 
  // corresponding c'tor:
  
  template < typename U ,
             typename V ,
             std::size_t Wt ,
             template < typename , std::size_t ,
                        typename , std::size_t > class W >
  simd_type ( const W < U , vsize , V , Wt > & ini )
  {
    *this = ini ;
  }
 
  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 ) ) ;
 
    static const index_type ix ( index_type::iota() ) ;
    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 (*this).
 
  // load uses a simple loop, which is about as easy to recognize as
  // an autovectorizable construct as it gets:
 
  void load ( const value_type * const p_src )
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      (*this) [ i ] = p_src [ i ] ;
  }
 
  // generic gather performs the gather operation using a loop.
  // Rather than loading consecutive values, The offset from 'p_src'
  // is looked up in 'indexes' for each vector element. We allow
  // any old indexable type as index_type, not just 'index_type'
  // defined above.
 
  template < typename index_type >
  void gather ( const value_type * const p_src ,
                const index_type & indexes )
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      (*this) [ i ] = p_src [ indexes [ i ] ] ;
  }
 
  // 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 container to memory
 
  void store ( value_type * const p_trg ) const
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      p_trg [ i ] = (*this) [ i ] ;
  }
 
  // scatter is the reverse operation to gather, see the comments there.
 
  template < typename index_type >
  void scatter ( value_type * const p_trg ,
                 const index_type & indexes ) const
  {
    for ( size_type i = 0 ; i < vsize ; i++ )
      p_trg [ indexes [ i ] ] = (*this) [ i ] ;
  }
 
  // '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 int & step )
  {
    index_type indexes ( IndexesFrom ( 0 , step ) ) ;
    gather ( p_src , indexes ) ;
  }
 
  void rscatter ( value_type * p_trg ,
                  const int & step ) const
  {
    index_type indexes ( IndexesFrom ( 0 , step ) ) ;
    scatter ( p_trg , indexes ) ;
  }
 
  // use 'indexes' to perform a gather from the data held in '(*this)'
  // and return the result of the gather operation.
 
  template < typename index_type >
  simd_type shuffle ( index_type indexes )
  {
    simd_type result ;
    for ( size_type i = 0 ; i < vsize ; i++ )
      result [ i ] = (*this) [ indexes [ i ] ] ;
    return result ;
  }
 
  // operator[] with an index_type argument performs the same
  // operation
 
  simd_type operator[] ( index_type indexes )
  {
    return shuffle ( 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 ) \
    { \
      simd_type result ; \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        result [ i ] = std::FUNC ( arg [ i ] ) ; \
      return result ; \
    }
 
  BROADCAST_STD_FUNC(abs)
  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)
 
  BROADCAST_STD_FUNC(sin)
  BROADCAST_STD_FUNC(cos)
  BROADCAST_STD_FUNC(tan)
  BROADCAST_STD_FUNC(asin)
  BROADCAST_STD_FUNC(acos)
  BROADCAST_STD_FUNC(atan)
 
  #undef BROADCAST_STD_FUNC
 
  #define BROADCAST_STD_FUNC2(FUNC) \
    friend simd_type FUNC ( simd_type arg1 , \
                            simd_type arg2 ) \
    { \
      simd_type result ; \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        result [ i ] = std::FUNC ( arg1 [ i ] , arg2 [ i ] ) ; \
      return result ; \
    }
 
  // 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 ) \
    { \
      simd_type result ; \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        result [ i ] = FUNC ( arg1 [ i ] , arg2 [ i ] , arg3[i] ) ; \
      return result ; \
    }
 
  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'
  // TODO: some operators make sense as simply manipulating the 'raw'
  // bits in the data (e.g. ~ or ^) - I currently limit them to int,
  // but this might be relaxed, casting the data to some unsigned
  // integral format of equal size.
 
  #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 itself are taken as rhs arguments.
 
  #define OPEQ_FUNC(OPFUNC,OPEQ,CONSTRAINT) \
    simd_type & OPFUNC ( value_type rhs ) \
    { \
      CONSTRAINT \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        (*this) [ i ] OPEQ rhs ; \
      return *this ; \
    } \
    simd_type & OPFUNC ( simd_type rhs ) \
    { \
      CONSTRAINT \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        (*this) [ i ] OPEQ rhs [ i ] ; \
      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 \
                       < is_scalar < RHST > :: value \
                       > :: type \
           > \
  simd_type < C_PROMOTE ( value_type , RHST ) , vsize > \
  OPFUNC ( simd_type < RHST , vsize > rhs ) const \
  { \
    CONSTRAINT \
    simd_type < C_PROMOTE ( value_type , RHST ) , vsize > help ( *this ) ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      help [ i ] = ( (*this) [ i ] OP rhs [ i ] ) ; \
    return help ; \
  } \
  template < typename RHST , \
             typename = typename std::enable_if \
                       < is_scalar < RHST > :: value \
                       > :: type \
           > \
  simd_type < C_PROMOTE ( value_type , RHST ) , vsize > \
  OPFUNC ( RHST rhs ) const \
  { \
    CONSTRAINT \
    simd_type < C_PROMOTE ( value_type , RHST ) , vsize > help ( *this ) ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      help [ i ] = ( (*this) [ i ] OP rhs ) ; \
    return help ; \
  } \
  template < typename LHST , \
             typename = typename std::enable_if \
                       < is_scalar < LHST > :: value \
                       > :: type \
           > \
  friend simd_type < C_PROMOTE ( LHST , value_type ) , vsize > \
  OPFUNC ( LHST lhs , simd_type rhs ) \
  { \
    CONSTRAINT \
    simd_type < C_PROMOTE ( LHST , value_type ) , vsize > help ( lhs ) ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      help [ i ] = ( lhs OP rhs [ i ] ) ; \
    return help ; \
  }
 
  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 \
    { \
      simd_type help ; \
      for ( size_type i = 0 ; i < vsize ; i++ ) \
        help [ i ] = OP (*this) [ i ] ; \
      return help ; \
    }
 
  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 ) \
  // { \
  //   mask_type m ; \
  //   for ( size_type i = 0 ; i < vsize ; i++ ) \
  //     m [ i ] = ( lhs [ i ] OP rhs [ i ] ) ; \
  //   return m ; \
  // } \
  // friend mask_type OPFUNC ( simd_type lhs , \
  //                           value_type rhs ) \
  // { \
  //   mask_type m ; \
  //   for ( size_type i = 0 ; i < vsize ; i++ ) \
  //     m [ i ] = ( lhs [ i ] OP rhs ) ; \
  //   return m ; \
  // } \
  // friend mask_type OPFUNC ( value_type lhs , \
  //                           simd_type rhs ) \
  // { \
  //   mask_type m ; \
  //   for ( size_type i = 0 ; i < vsize ; i++ ) \
  //     m [ i ] = ( lhs OP rhs [ i ] ) ; \
  //   return m ; \
  // }
 
#define COMPARE_FUNC(OPFUNC,OP) \
  template < typename RHST , \
             typename = typename std::enable_if \
                       < is_scalar < RHST > :: value \
                       > :: type \
           > \
  mask_type OPFUNC ( simd_type < RHST , vsize > rhs ) const \
  { \
    mask_type m ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      m [ i ] = ( (*this) [ i ] OP rhs [ i ] ) ; \
    return m ; \
  } \
  template < typename RHST , \
             typename = typename std::enable_if \
                       < is_scalar < RHST > :: value \
                       > :: type \
           > \
  mask_type OPFUNC ( RHST rhs ) const \
  { \
    mask_type m ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      m [ i ] = ( (*this) [ i ] OP rhs ) ; \
    return m ; \
  } \
  template < typename LHST , \
             typename = typename std::enable_if \
                       < is_scalar < LHST > :: value \
                       > :: type \
           > \
  friend mask_type OPFUNC ( LHST lhs , simd_type rhs ) \
  { \
    mask_type m ; \
    for ( size_type i = 0 ; i < vsize ; i++ ) \
      m [ i ] = ( lhs OP rhs [ i ] ) ; \
    return m ; \
  }
 
  COMPARE_FUNC(operator<,<) ;
  COMPARE_FUNC(operator<=,<=) ;
  COMPARE_FUNC(operator>,>) ;
  COMPARE_FUNC(operator>=,>=) ;
  COMPARE_FUNC(operator==,==) ;
  COMPARE_FUNC(operator!=,!=) ;
 
  #undef COMPARE_FUNC
 
  // next we define a masked vector as an object holding two pieces of
  // information: one mask, determining which of the vector's elements
  // will be 'open' to an effect, and one reference to a vector, which
  // will be affected by an assignment or augmented assignment to the
  // masked_type object. The mask is not held by reference - it is
  // typically created right inside an invocation of the V(M) = ...
  // idiom and will only persist when copied to the masked_type object.
  // The compiler will take care of the masks's lifetime and storage and
  // make sure the process is efficient and avoids unnecessary copying.
  // The resulting object will only be viable as long as the referred-to
  // vector is kept '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 - the V(M)=...
  // idiom mentioned above.
  // Note that this does not have any effect on those values in 'whither'
  // for which the mask is false. They remain unchanged.
 
  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 \
        for ( size_type i = 0 ; i < vsize ; i++ ) \
        { \
          if ( whether [ i ] ) \
            whither [ i ] OPEQ rhs ; \
        } \
        return whither ; \
      } \
      simd_type & OPFUNC ( simd_type rhs ) \
      { \
        CONSTRAINT \
        for ( size_type i = 0 ; i < vsize ; i++ ) \
        { \
          if ( whether [ i ] ) \
            whither [ i ] OPEQ rhs [ i ] ; \
        } \
        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 ) ;
  }
 
  // next we have a few functions creating masks - mainly for Vc
  // compatibility.
 
  static mask_type isnegative ( const simd_type & rhs )
  {
    return ( rhs < value_type(0) ) ;
  }
 
  static mask_type isfinite ( const simd_type & rhs )
  {
    for ( std::size_t i = 0 ; i < vsize ; i++ )
    {
      if ( isInf ( rhs[i] ) )
        return false ;
    }
    return true ;
  }
 
  static mask_type isnan ( const simd_type & rhs )
  {
    for ( std::size_t i = 0 ; i < vsize ; i++ )
    {
      if ( isNan ( rhs[i] ) )
        return true ;
    }
    return false ;
  }
 
  // 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 \
    { \
      simd_type result ( threshold ) ; \
      for ( std::size_t i = 0 ; i < vsize ; i++ ) \
      { \
        if ( (*this) [ i ] REL threshold ) \
          result [ i ] = (*this) [ i ] ; \
      } \
      return result ; \
    }
 
  CLAMP(at_least,>)
  CLAMP(at_most,<)
 
  #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 ;
  }
} ;
 
// reductions for masks. It's often necessary to determine whether
// a mask is completely full or empty, or has at least some non-false
// members. We allow arbitrary simd_type to be testes with any_of,
// all_of and none_of, rather than enforcing that only 'proper' masks
// (namely simd_type < bool , vsize >) can be tested.
 
template < typename T , std::size_t vsize >
bool any_of ( simd_type < T , vsize > arg )
{
  bool result = false ;
  for ( std::size_t i = 0 ; i < vsize ; i++ )
    result |= bool ( arg [ i ] ) ;
  return result ;
}
 
template < typename T , std::size_t vsize >
bool all_of ( simd_type < T , vsize > arg )
{
  bool result = true ;
  for ( std::size_t i = 0 ; i < vsize ; i++ )
    result &= bool ( arg [ i ] ) ;
  return result ;
}
 
template < typename T , std::size_t vsize >
bool none_of ( simd_type < T , vsize > arg )
{
  bool result = true ;
  for ( std::size_t i = 0 ; i < vsize ; i++ )
    result &= ( ! bool ( arg [ i ] ) ) ;
  return result ;
}
 
} ;
 
template < typename T , std::size_t N >
struct std::allocator_traits < vspline::simd_type < T , N > >
{
  typedef std::allocator < vspline::simd_type < T , N > > type ;
} ;
 
#endif // #define VSPLINE_SIMD_TYPE_H