text: /************************************************************************/
  /*                                                                      */
  /*               Copyright 2004-2005 by Ullrich Koethe                  */
  /*                                                                      */
  /*    This file is part of the VIGRA computer vision library.           */
  /*    The VIGRA Website is                                              */
  /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
  /*    Please direct questions, bug reports, and contributions to        */
  /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
  /*        vigra@informatik.uni-hamburg.de                               */
  /*                                                                      */
  /*    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.                                   */                
  /*                                                                      */
  /************************************************************************/
  
  #ifndef VIGRA_FIXEDPOINT_HXX
  #define VIGRA_FIXEDPOINT_HXX
  
  #include "mathutil.hxx"
  #include "static_assert.hxx"
  #include "error.hxx"
  #include "numerictraits.hxx"
  
  namespace vigra {
  
  template <unsigned IntBits, unsigned FractionalBits>
  class FixedPoint;
  
  struct Error_FixedPointTraits_not_specialized_for_this_case;
  
  template <class T1, class T2>
  class FixedPointTraits
  {
  public:
      typedef Error_FixedPointTraits_not_specialized_for_this_case PlusType;
      typedef Error_FixedPointTraits_not_specialized_for_this_case MinusType;
      typedef Error_FixedPointTraits_not_specialized_for_this_case MultipliesType;
  //    typedef Error_FixedPointTraits_not_specialized_for_this_case DividesType;
  };
  
  // return type policy: 
  //     * try to allocate enough bits to represent the biggest possible result
  //     * in case of add/subtract: if all bits of the internal int are used up, 
  //                                keep the representation
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
  {
      enum { MaxIntBits  = (IntBits1 < IntBits2) ? IntBits2 : IntBits1,
             MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1,
             PlusMinusIntBits = (MaxIntBits + 1 + MaxFracBits < 32) ?
                                 MaxIntBits + 1 : MaxIntBits,
             MultipliesFracBits = (IntBits1 + IntBits2 < 31) 
                                      ? (FracBits1 + FracBits2) > (31 - IntBits1 - IntBits2) 
                                              ? 31 - IntBits1 - IntBits2
                                              : FracBits1 + FracBits2
                                      : 0
           };
  public:
      typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
      typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
      typedef FixedPoint<IntBits1 + IntBits2, MultipliesFracBits>  MultipliesType;
  //    typedef FixedPoint<IntBits1 + FracBits2, FracBits1 + IntBits2>  DividesType;
  };
  
  template <unsigned IntBits, unsigned FracBits>
  struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
  {
      enum { SRTotalBits = (IntBits + FracBits + 1) / 2,
             SRIntBits   = (IntBits + 1) / 2,
             SRFracBits  = SRTotalBits - SRIntBits
           };
  public:
      typedef FixedPoint<IntBits, FracBits>      Type;
      typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
      typedef Type                               SquareRootArgument;
  };
  
  
  #ifndef DOXYGEN
  
  template <int N>
  struct FixedPoint_overflow_error__More_than_31_bits_requested
  : staticAssert::AssertBool<(N < 32)>
  {};
  
  #endif /* DOXYGEN */
  
  
  
  template <bool Predicate>
  struct FixedPoint_assignment_error__Target_object_has_too_few_integer_bits
  : staticAssert::AssertBool<Predicate>
  {};
  
  enum FixedPointNoShift { FPNoShift };
  
  namespace detail {
  
  template <bool MustRound>
  struct FPAssignWithRound;
  
  template <>
  struct FPAssignWithRound<false>
  {
      template <int N>
      static inline int exec(int v) { return v << (-N); }
  };
  
  template <>
  struct FPAssignWithRound<true>
  {
      template <int N>
      static inline int exec(int const v)
      {
          return (v + (1 << (N - 1))) >> (N);
      }
  };
  
  template <bool MustRound>
  struct FPMulImplementation;
  
  template <>
  struct FPMulImplementation<false>
  {
      template <int N>
      static inline int exec(int l, int r) { return (l * r) << (-N); }
  };
  
  template <>
  struct FPMulImplementation<true>
  {
      template <int N>
      static inline int exec(int l, int r)
      {
          // there is not enough space in the result
          // => perform calculations that preserve as much accuracy as possible
          enum { diffl = N / 2, diffr = N  - diffl, maskl = (1 << diffl) - 1, maskr = (1 << diffr) - 1 };
          int shiftl = l >> diffl;
          int shiftr = r >> diffr;
  
          return shiftl * shiftr + (((l & maskl) * shiftr) >> diffl) +
                                   (((r & maskr) * shiftl) >> diffr);
      }
  };
  
  } // namespace detail
  
  /********************************************************/
  /*                                                      */
  /*                      FixedPoint                      */
  /*                                                      */
  /********************************************************/
  
  /** Template for fixed point arithmetic.
  
      Fixed point arithmetic is used when computations with fractional accuracy
      must be made at the highest speed possible (e.g. in the inner loop
      of a volume rendering routine). The speed-up relative to floating
      point arithmetic can be dramatic, especially when one can avoid
      conversions between integer and floating point numbers (these are 
      very expensive because integer and floating point arithmetic
      resides in different pipelines). 
      
      The template wraps an <tt>int</tt> and uses <tt>IntBits</tt> to
      represent the integral part of a number, and <tt>FractionalBits</tt>
      for the fractional part, where <tt>IntBits + FractionalBits < 32</tt>.
      (The 32rd bit is reserved because FixedPoint is a signed type).
      These numbers will be automatically allocated in an intelligent way
      in the result of an arithmetic operation. For example, when two 
      fixed point numbers are multiplied, the required number of integer
      bits in the result is the sum of the number of integer bits of the
      arguments, but only when so many bits are available. This is figured out
      by means of FixedPointTraits, and a compile-time error is raised
      when no suitable representation can be found. The idea is that the right
      thing happens automatically as often as possible.
  
      <tt>FixedPoint</tt> implements the required interface of an
      \ref AlgebraicRing and the required numeric and
      promotion traits. In addition, it supports functions <tt>add</tt>, 
      <tt>sub</tt>, and <tt>mul</tt>, where a particular layout of the result can
      be enforced. 
      
      <tt>unsigned char, signed char, unsigned short, signed short, int</tt> can be
      transformed into a FixedPoint with appropriate layout by means of the factory
      function <tt>fixedPoint()</tt>.
  
      <b>See also:</b>
      <ul>
      <li> \ref FixedPointOperations
      <li> \ref FixedPointTraits
      </ul>
  
      <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
      Namespace: vigra
  */
  template <unsigned IntBits, unsigned FractionalBits>
  class FixedPoint
  {
  public:
      enum {
          INT_BITS        = IntBits,
          FRACTIONAL_BITS = FractionalBits,
          TOTAL_BITS      = IntBits + FractionalBits,
          MAX             = (int)(((unsigned)1 << TOTAL_BITS) - 1),
          ONE             = 1 << FractionalBits,
          ONE_HALF        = ONE >> 1,
          FRACTIONAL_MASK = ONE - 1,
          INT_MASK        = MAX ^ FRACTIONAL_MASK
      };
  
      Int32 value;
  
      FixedPoint()
      {
          VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
      }
  
          /** Construct from an int (fractional part will become zero).
          */
      explicit FixedPoint(int v)
      : value(v << FractionalBits)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
      }
  
          /** Construct from an int by a bitwise copy. This is normally only used internally.
          */
      FixedPoint(int v, FixedPointNoShift)
      : value(v)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
      }
  
          /** Construct from an double and round the fractional part to 
              <tt>FractionalBits</tt> accuracy. A PreconditionViolation exception is raised when
              the integer part is too small to represent the number.
          */
      explicit FixedPoint(double rhs)
      : value((int)round(rhs * ONE))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
          vigra_precondition(abs(rhs * ONE) <= (double)MAX,
              "FixedPoint(double rhs): Too few integer bits to convert rhs.");
      }
  
          /** Copy constructor.
          */
      FixedPoint(const FixedPoint &other)
      : value(other.value)
      {}
  
          /** Construct from a FixedPoint with different layout. It rounds as appropriate and raises
              a compile-time error when the target type has too few integer bits.
          */
      template <unsigned Int2, unsigned Frac2>
      FixedPoint(const FixedPoint<Int2, Frac2> &other)
      : value(detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
          VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
      }
  
          /** Assignment from int. The fractional part will become zero.  
              A PreconditionViolation exception is raised when
              the integer part is too small to represent the number.
          */
      FixedPoint &operator=(int rhs)
      {
          vigra_precondition(abs(rhs) < (1 << IntBits),
              "FixedPoint::operator=(int rhs): Too few integer bits to represent rhs.");
          value = rhs << FractionalBits;
          return *this;
      }
  
          /** Assignment form double. The fractional part is rounded, and a 
              PreconditionViolation exception is raised when
              the integer part is too small to represent the number.
          */
      FixedPoint &operator=(double rhs)
      {
          vigra_precondition(abs(rhs) <= ((1 << IntBits) - 1),
              "FixedPoint::operator=(double rhs): Too few integer bits to convert rhs.");
          value = (int)round(rhs * ONE);
          return *this;
      }
  
          /** Copy assignment.
          */
      FixedPoint & operator=(const FixedPoint &other)
      {
          value = other.value;
          return *this;
      }
  
          /** Assignment from a FixedPoint with different layout. It rounds as appropriate and raises
              a compile-time error when the target type has too few integer bits.
          */
      template <unsigned Int2, unsigned Frac2>
      FixedPoint & operator=(const FixedPoint<Int2, Frac2> &other)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
          value = detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value);
          return *this;
      }
  
          /** Negation.
          */
      FixedPoint operator-() const
      {
          return FixedPoint(-value, FPNoShift);
      }
  
          /** Pre-increment.
          */
      FixedPoint & operator++()
      {
          value += ONE;
          return *this;
      }
  
          /** Post-increment.
          */
      FixedPoint operator++(int)
      {
          FixedPoint old(*this);
          value += ONE;
          return old;
      }
  
          /** Pre-decrement.
          */
      FixedPoint & operator--()
      {
          value -= ONE;
          return *this;
      }
  
          /** Post-decrement.
          */
      FixedPoint operator--(int)
      {
          FixedPoint old(*this);
          value -= ONE;
          return old;
      }
  
          /** Add-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
              a compile-time error when the target type has too few integer bits.
          */
      template <unsigned Int2, unsigned Frac2>
      FixedPoint & operator+=(const FixedPoint<Int2, Frac2> &other)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
          value += detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
          return *this;
      }
  
          /** Subtract-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
              a compile-time error when the target type has too few integer bits.
          */
      template <unsigned Int2, unsigned Frac2>
      FixedPoint & operator-=(const FixedPoint<Int2, Frac2> &other)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
          value -= detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
          return *this;
      }
      
          /** Multiply-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
              a compile-time error when the target type has too few integer bits.
          */
      template <unsigned Int2, unsigned Frac2>
      FixedPoint & operator*=(const FixedPoint<Int2, Frac2> &other)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
          value = detail::FPMulImplementation<(Frac2 > 0)>::template exec<Frac2>(value, other.value);
          return *this;
      }
  };
  
  #define VIGRA_FIXED_POINT_FACTORY(T, INTBITS) \
      inline FixedPoint<INTBITS, 0> fixedPoint(T t) \
      { \
          return FixedPoint<INTBITS, 0>(t, FPNoShift); \
      }
  
  VIGRA_FIXED_POINT_FACTORY(unsigned char, 8)
  VIGRA_FIXED_POINT_FACTORY(signed char, 7)
  VIGRA_FIXED_POINT_FACTORY(unsigned short, 16)
  VIGRA_FIXED_POINT_FACTORY(signed short, 15)
  VIGRA_FIXED_POINT_FACTORY(int, 31)
  
  #undef VIGRA_FIXED_POINT_FACTORY
  
  template <class T>
  struct FixedPointCast;
  
  #define VIGRA_FIXED_POINT_CAST(type) \
  template <> \
  struct FixedPointCast<type> \
  { \
      template <unsigned IntBits, unsigned FracBits> \
      static type cast(FixedPoint<IntBits, FracBits> v) \
      { \
          return round(v); \
      } \
  };
  
  VIGRA_FIXED_POINT_CAST(Int8)
  VIGRA_FIXED_POINT_CAST(UInt8)
  VIGRA_FIXED_POINT_CAST(Int16)
  VIGRA_FIXED_POINT_CAST(UInt16)
  VIGRA_FIXED_POINT_CAST(Int32)
  VIGRA_FIXED_POINT_CAST(UInt32)
  
  #undef VIGRA_FIXED_POINT_CAST
  
  template <>
  struct FixedPointCast<float>
  {
      template <unsigned IntBits, unsigned FracBits>
      static float cast(FixedPoint<IntBits, FracBits> v)
      {
          return (float)v.value / FixedPoint<IntBits, FracBits>::ONE;
      }
  };
  
  template <>
  struct FixedPointCast<double>
  {
      template <unsigned IntBits, unsigned FracBits>
      static double cast(FixedPoint<IntBits, FracBits> v)
      {
          return (double)v.value / FixedPoint<IntBits, FracBits>::ONE;
      }
  };
  
  /********************************************************/
  /*                                                      */
  /*                 FixedPointOperations                 */
  /*                                                      */
  /********************************************************/
  
  /** \addtogroup FixedPointOperations Functions for FixedPoint
  
      \brief     <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
  
      These functions fulfill the requirements of an \ref AlgebraicRing.
  
      Namespace: vigra
      <p>
  
   */
  //@{
  
      /** Convert a FixedPoint to a built-in type.
          If the target is integral, the value is rounded.<br>
          Usage:
          \code
          FixedPoint<16,15> fp(...);
          
          double d = fixed_point_cast<double>(fp);
          \endcode
      */
  template <class TARGET, unsigned IntBits, unsigned FracBits>
  TARGET fixed_point_cast(FixedPoint<IntBits, FracBits> v)
  {
      return FixedPointCast<TARGET>::cast(v);
  }
  
      /// equal
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator==(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) == (r.value << (MaxFracBits - FracBits2));
  }
  
      /// not equal
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator!=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) != (r.value << (MaxFracBits - FracBits2));
  }
  
      /// less than
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator<(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) < (r.value << (MaxFracBits - FracBits2));
  }
  
      /// less or equal
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator<=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) <= (r.value << (MaxFracBits - FracBits2));
  }
  
      /// greater
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator>(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) > (r.value << (MaxFracBits - FracBits2));
  }
  
      /// greater or equal
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  bool operator>=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return (l.value << (MaxFracBits - FracBits1)) >= (r.value << (MaxFracBits - FracBits2));
  }
  
      /// addition with automatic determination of the appropriate result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
  operator+(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return typename
          FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
          PlusType((l.value << (MaxFracBits - FracBits1)) + (r.value << (MaxFracBits - FracBits2)), FPNoShift);
  }
  
      /// addition with enforced result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
            unsigned IntBits3, unsigned FracBits3>
  inline void
  add(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
      FixedPoint<IntBits3, FracBits3> & result)
  {
      result = l + r;
  }
  
      /// subtraction with automatic determination of the appropriate result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MinusType
  operator-(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
      return typename
          FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
          MinusType((l.value << (MaxFracBits - FracBits1)) - (r.value << (MaxFracBits - FracBits2)), FPNoShift);
  }
  
      /// subtraction with enforced result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
            unsigned IntBits3, unsigned FracBits3>
  inline void
  sub(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
      FixedPoint<IntBits3, FracBits3> & result)
  {
      result = l - r;
  }
  
      /// multiplication with automatic determination of the appropriate result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  inline
  typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MultipliesType
  operator*(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
  {
      typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
          MultipliesType res;
      mul(l, r, res);
      return res;
  }
  
      /// multiplication with enforced result type.
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
            unsigned IntBits3, unsigned FracBits3>
  inline void
  mul(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
      FixedPoint<IntBits3, FracBits3> & result)
  {
      VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits1 + IntBits2 <= IntBits3)>));
      enum { diff = FracBits1 + FracBits2 - FracBits3 };
      result.value = detail::FPMulImplementation<(diff > 0)>::template exec<diff>(l.value, r.value);
  }
  
      /// square root.
  template <unsigned IntBits, unsigned FracBits>
  inline typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult
  sqrt(FixedPoint<IntBits, FracBits> v)
  {
      return typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult(sqrti(v.value), FPNoShift);
  }
  
      /// absolute value.
  template <unsigned IntBits, unsigned FracBits>
  inline FixedPoint<IntBits, FracBits>
  abs(FixedPoint<IntBits, FracBits> v)
  {
      return FixedPoint<IntBits, FracBits>(abs(v.value), FPNoShift);
  }
  
      /// squared norm (same as v*v).
  template <unsigned IntBits, unsigned FracBits>
  inline
  typename FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
  squaredNorm(FixedPoint<IntBits, FracBits> v)
  {
      return v*v;
  }
  
      /// norm (same as abs).
  template <unsigned IntBits, unsigned FracBits>
  inline
  FixedPoint<IntBits, FracBits>
  norm(FixedPoint<IntBits, FracBits> const & v)
  {
      return abs(v);
  }
  
      /// fractional part.
  template <unsigned IntBits, unsigned FracBits>
  inline FixedPoint<0, FracBits>
  frac(FixedPoint<IntBits, FracBits> v)
  {
      return FixedPoint<0, FracBits>(v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK, FPNoShift);
  }
  
      /// dual fractional part: <tt>1 - frac(v)</tt>.
  template <unsigned IntBits, unsigned FracBits>
  inline FixedPoint<0, FracBits>
  dual_frac(FixedPoint<IntBits, FracBits> v)
  {
      return FixedPoint<0, FracBits>(FixedPoint<0, FracBits>::ONE - 
                                     (v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK), FPNoShift);
  }
  
      /// rounding down.
  template <unsigned IntBits, unsigned FracBits>
  inline int
  floor(FixedPoint<IntBits, FracBits> v)
  {
      return(v.value >> FracBits);
  }
  
      /// rounding up.
  template <unsigned IntBits, unsigned FracBits>
  inline int
  ceil(FixedPoint<IntBits, FracBits> v)
  {
      return((v.value + FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK) >> FracBits);
  }
  
      /// rounding to the nearest integer.
  template <unsigned IntBits, unsigned FracBits>
  inline int
  round(FixedPoint<IntBits, FracBits> v)
  {
      return((v.value + FixedPoint<IntBits, FracBits>::ONE_HALF) >> FracBits);
  }
  
  //@}
  
  /********************************************************/
  /*                                                      */
  /*                     FixedPoint-Traits                */
  /*                                                      */
  /********************************************************/
  
  /** \page FixedPointTraits Numeric and Promote Traits of FixedPoint
  
      The numeric and promote traits for FixedPoint follow
      the general specifications for \ref NumericPromotionTraits and
      \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
      partial template specialization:
  
      \code
  
      template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
      class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
      {
          typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
          typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
          typedef FixedPoint<IntBits1 + IntBits2, FracBits1 + FracBits2>  MultipliesType;
      };
  
      template <unsigned IntBits, unsigned FracBits>
      struct NumericTraits<FixedPoint<IntBits, FracBits> >
      {
          typedef FixedPoint<IntBits, FracBits> Type;
              // Promote undefined because it depends on the layout, use FixedPointTraits
              // RealPromote in AlgebraicRing -- multiplication with double is not supported.
              // ComplexPromote in AlgebraicRing -- multiplication with double is not supported.
          typedef Type ValueType;
  
          typedef VigraFalseType isIntegral;
          typedef VigraTrueType  isScalar;
          typedef VigraTrueType  isSigned;
          typedef VigraTrueType  isOrdered;
          typedef VigraFalseType isComplex;
  
          ... // etc.
      };
  
      template <unsigned IntBits, unsigned FracBits>
      struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
      {
          typedef FixedPoint<IntBits, FracBits>      Type;
          typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
          typedef Type                               SquareRootArgument;
      };
      
      template <unsigned IntBits, unsigned FracBits>
      struct NormTraits<FixedPoint<IntBits, FracBits> >
      {
          typedef FixedPoint<IntBits, FracBits>         Type;
          typedef typename 
              FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
                                                        SquaredNormType;
          typedef Type                                  NormType;
      };
  
      template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
      struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
                           FixedPoint<IntBits2, FracBits2> >
      {
          typedef typename 
              FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType 
              Promote;
      };
      \endcode
  
      <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
      Namespace: vigra
  
  */
  template <unsigned IntBits, unsigned FracBits>
  struct NumericTraits<FixedPoint<IntBits, FracBits> >
  {
      typedef FixedPoint<IntBits, FracBits> Type;
          //typedef FixedPoint<IntBits, FracBits> Promote;
          //typedef FixedPoint<IntBits, FracBits> RealPromote;
          //typedef std::complex<RealPromote> ComplexPromote;
      typedef Type ValueType;
  
      typedef VigraFalseType isIntegral;
      typedef VigraTrueType  isScalar;
      typedef VigraTrueType  isSigned;
      typedef VigraTrueType  isOrdered;
      typedef VigraFalseType isComplex;
  
      static Type zero() { return Type(0, FPNoShift); }
      static Type one() { return Type(Type::ONE, FPNoShift); }
      static Type nonZero() { return one(); }
      static Type epsilon() { return Type(1, FPNoShift); }
      static Type smallestPositive() { return Type(1, FPNoShift); }
      static Type max() { return Type( Type::MAX, FPNoShift); }
      static Type min() { return -max(); }
  };
  
  template <unsigned IntBits, unsigned FracBits>
  struct NormTraits<FixedPoint<IntBits, FracBits> >
  {
      typedef FixedPoint<IntBits, FracBits>         Type;
      typedef typename 
          FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
                                                    SquaredNormType;
      typedef Type                                  NormType;
  };
  
  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
  struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
                       FixedPoint<IntBits2, FracBits2> >
  {
      typedef typename 
          FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType 
          Promote;
  };
  
  /***********************************************************************************/
  
  enum FPOverflowHandling { FPOverflowIgnore, FPOverflowSaturate, FPOverflowError };
  
  template <int IntBits, FPOverflowHandling OverflowHandling = FPOverflowIgnore>
  class FixedPoint16;
  
  /********************************************************/
  /*                                                      */
  /*                     FixedPoint16-Traits              */
  /*                                                      */
  /********************************************************/
  
  /** \page FixedPoint16Traits Numeric and Promote Traits of FixedPoint16
  
      The numeric and promote traits for FixedPoint16 follow
      the general specifications for \ref NumericPromotionTraits and
      \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
      partial template specialization:
  
      \code
      template <int IntBits, FPOverflowHandling OverflowHandling>
      struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
      {
          typedef FixedPoint16<IntBits, OverflowHandling> Type;
          typedef Type                                    Promote;
              // RealPromote undefined -- multiplication with double is not supported.
              // ComplexPromote undefined -- multiplication with double is not supported.
          typedef Type ValueType;
  
          typedef VigraFalseType isIntegral;
          typedef VigraTrueType  isScalar;
          typedef VigraTrueType  isSigned;
          typedef VigraTrueType  isOrdered;
          typedef VigraFalseType isComplex;
  
          ... // etc.
      };
  
      template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
      struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
                           FixedPoint16<IntBits2, OverflowHandling> >
      {
          typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
          ... // etc.
      };
  
      template <int IntBits, FPOverflowHandling OverflowHandling>
      struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
      {
          typedef FixedPoint16<IntBits, OverflowHandling>     Type;
          typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
          typedef Type                                        NormType;
      };
  
      template <int IntBits, FPOverflowHandling OverflowHandling>
      struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
      {
          typedef FixedPoint16<IntBits, OverflowHandling>            Type;
          typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
          typedef Type                                               SquareRootArgument;
      };
      \endcode
  
      <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
      Namespace: vigra
  
  */
  template <int IntBits, FPOverflowHandling OverflowHandling>
  struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
  {
      typedef FixedPoint16<IntBits, OverflowHandling> Type;
      typedef Type                                    Promote;
          // RealPromote undefined -- multiplication with double is not supported.
          // ComplexPromote undefined -- multiplication with double is not supported.
      typedef Type ValueType;
  
      typedef VigraFalseType isIntegral;
      typedef VigraTrueType  isScalar;
      typedef VigraTrueType  isSigned;
      typedef VigraTrueType  isOrdered;
      typedef VigraFalseType isComplex;
  
      static Type zero() { return Type(0, FPNoShift); }
      static Type one() { return Type(Type::ONE, FPNoShift); }
      static Type nonZero() { return one(); }
      static Type epsilon() { return Type(1, FPNoShift); }
      static Type smallestPositive() { return Type(1, FPNoShift); }
      static Type max() { return Type( Type::MAX, FPNoShift); }
      static Type min() { return Type( Type::MIN, FPNoShift); }
  
      static Promote toPromote(Type v) { return v; }
      static Type fromPromote(Promote v) { return v; }; 
  };
  
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
                       FixedPoint16<IntBits2, OverflowHandling> >
  {
      typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
      static Promote toPromote(FixedPoint16<IntBits1, OverflowHandling> v) { return Promote(v); }
      static Promote toPromote(FixedPoint16<IntBits2, OverflowHandling> v) { return Promote(v); }
  };
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  struct PromoteTraits<FixedPoint16<IntBits, OverflowHandling>,
                       FixedPoint16<IntBits, OverflowHandling> >
  {
      typedef FixedPoint16<IntBits, OverflowHandling> Promote;
      static Promote toPromote(FixedPoint16<IntBits, OverflowHandling> v) { return v; }
  };
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
  {
      typedef FixedPoint16<IntBits, OverflowHandling>     Type;
      typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
      typedef Type                                        NormType;
  };
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
  {
      typedef FixedPoint16<IntBits, OverflowHandling>            Type;
      typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
      typedef Type                                               SquareRootArgument;
  };
  
  #ifndef DOXYGEN
  
  template <bool Compatible>
  struct FixedPoint_error__Right_shift_operator_has_unsupported_semantics
  : staticAssert::AssertBool<Compatible>
  {};
  
  #endif /* DOXYGEN */
  
  template <bool Predicate>
  struct FixedPoint16_assignment_error__Target_object_has_too_few_integer_bits
  : staticAssert::AssertBool<Predicate>
  {};
  
  namespace detail {
  
  template<int BeforeIntBits, int AfterIntBits, 
           bool Round = false,
           bool RightShift = (AfterIntBits >= BeforeIntBits)>
  struct FP16Align;
  
  template<int BeforeIntBits>
  struct FP16Align<BeforeIntBits, BeforeIntBits, true, true>
  {
      static inline Int32 exec(Int32 v)
      {
          return v;
      }
  };
  
  template<int BeforeIntBits>
  struct FP16Align<BeforeIntBits, BeforeIntBits, false, true>
  {
      static inline Int32 exec(Int32 v)
      {
          return v;
      }
  };
  
  template<int BeforeIntBits, int AfterIntBits>
  struct FP16Align<BeforeIntBits, AfterIntBits, false, true>
  {
      static inline Int32 exec(Int32 v)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
          return v >> (AfterIntBits - BeforeIntBits);
      }
  };
  
  template<int BeforeIntBits, int AfterIntBits>
  struct FP16Align<BeforeIntBits, AfterIntBits, true, true>
  {
      enum { ONE_HALF = 1 << (AfterIntBits - BeforeIntBits - 1) };
      static inline Int32 exec(Int32 v)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
          return (v + ONE_HALF) >> (AfterIntBits - BeforeIntBits);
      }
  };
  
  template<int BeforeIntBits, int AfterIntBits, bool Round>
  struct FP16Align<BeforeIntBits, AfterIntBits, Round, false>
  {
      static inline Int32 exec(Int32 v)
      {
          return v << (BeforeIntBits - AfterIntBits);
      }
  };
  
  template <FPOverflowHandling OverflowHandling = FPOverflowIgnore>
  struct FP16OverflowHandling
  {
      static inline Int32 exec(Int32 v)
      {
          return v;
      }
      
      static inline Int32 exec(UInt32 v)
      {
          return v;
      }
  };
  
  template <>
  struct FP16OverflowHandling<FPOverflowSaturate>
  {
      static inline Int32 exec(Int32 v)
      {
          if(v >= 1 << 15)
              return (1 << 15) - 1;
          if(v < -(1 << 15))
              return -(1 << 15);
          return v;
      }
      static inline Int32 exec(UInt32 v)
      {
          if(v >= 1 << 15)
              return (1 << 15) - 1;
          return v;
      }
  };
  
  template <>
  struct FP16OverflowHandling<FPOverflowError>
  {
      static inline Int32 exec(Int32 v)
      {
          vigra_precondition(v < (1 << 15) && v >= -(1 << 15),
               "FixedPoint16: Operation overflows.");
          return v;
      }
      static inline Int32 exec(UInt32 v)
      {
          vigra_precondition(v < (1 << 15),
               "FixedPoint16: Operation overflows.");
          return v;
      }
  };
  
  
  template <int IntBits1, int IntBits2, int IntBitsOut, 
            FPOverflowHandling OverflowHandling >
  struct FP16AddImpl
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      static inline Int32 exec(Int32 t1, Int32 t2)
      {
         return FP16OverflowHandling<OverflowHandling>::exec(
                    FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
                        FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) + 
                        FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
      }
  };
  
  template <int IntBits1, int IntBits2, int IntBitsOut, 
            FPOverflowHandling OverflowHandling >
  struct FP16SubImpl
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      static inline Int32 exec(Int32 t1, Int32 t2)
      {
         return FP16OverflowHandling<OverflowHandling>::exec(
                    FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
                        FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) - 
                        FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
      }
  };
  
  template <int IntBits1, int IntBits2, int IntBitsOut, 
            FPOverflowHandling OverflowHandling >
  struct FP16MulImpl
  {
      static inline Int32 exec(Int32 t1, Int32 t2)
      {
          return FP16OverflowHandling<OverflowHandling>::exec(
                     FP16Align<IntBits1+IntBits2, IntBitsOut+15, /*Round*/ true>::exec(t1*t2));
      }
  };
  
  template <int IntBits1, int IntBits2, int IntBitsOut, 
            FPOverflowHandling OverflowHandling >
  struct FP16DivImpl
  {
      static inline Int32 exec(Int32 t1, Int32 t2)
      {
          if(t2 == 0)
              return (t1 >= 0)
                         ?  (1 << 15) - 1
                         : -(1 << 15);
          return FP16OverflowHandling<OverflowHandling>::exec(
                     FP16Align<IntBits1-IntBits2, IntBitsOut+1, /*Round*/ true>::exec((t1<<16)/t2));
      }
  };
  
  } // namespace detail
  
  /********************************************************/
  /*                                                      */
  /*                      FixedPoint16                    */
  /*                                                      */
  /********************************************************/
  
  template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
  TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v);
  
  /** Template for 16-bit signed fixed point arithmetic.
  
      Fixed point arithmetic is used when computations with fractional accuracy
      must be made at the highest speed possible (e.g. in the inner loop
      of a volume rendering routine). The speed-up relative to floating
      point arithmetic can be dramatic, especially when one can avoid
      conversions between integer and floating point numbers (these are 
      very expensive because integer and floating point arithmetic
      resides in different pipelines). 
      
      The template wraps an <tt>Int16</tt> and uses <tt>IntBits</tt> to
      represent the integral part of a number, and <tt>15 - IntBits</tt>
      for the fractional part. The 16th bit is reserved because FixedPoint16 
      is a signed type. Results of expressions with mixed types will preserve
      larger number of <tt>IntBits</tt> of the results, in order to minimize
      the possibility for overflow. Nonetheless, overflow can occur, and the 
      template parameter <tt>OverflowHandling</tt> determines how this will be
      handled:
      
      <DL>
      <DT>FPOverflowIgnore<DD> (default) Ignore overflow, i.e. use the usual modulo behavior of the
                               built-in integer types.
                         
      <DT>FPOverflowSaturate<DD> Use the largest or smallest representable number (depending on sign)
                                 in case of overflow.
  
      <DT>FPOverflowError<DD> Throw <tt>PreconditionViolation</tt> upon overflow. This is useful for 
                              debugging.
      </DL>
      
      The implementation relies on Int32-arithmetic and requires that the right-shift operator
      preserves signedness. Although not enforced by the C++ standard, this is implemented
      by most of today's processors. This property is checked by a 
      VIGRA_STATIC_ASSERT(FixedPoint_error__Right_shift_operator_has_unsupported_semantics).
  
      <tt>FixedPoint16</tt> implements the required interface of an
      \ref AlgebraicRing and the required numeric and
      promotion traits. In addition, it supports functions <tt>add</tt>, 
      <tt>sub</tt>, <tt>mul</tt>, and <tt>div</tt>, where a particular layout 
      of the result can be enforced. 
      
      Built-in numeric types can be converted into <tt>FixedPoint16</tt> by the 
      appropriate constructors, and from <tt>FixedPoint16</tt> by means of
      <tt>fixed_point_cast&lt;TargetType&gt;(fixedPoint)</tt>.
  
      <b>See also:</b>
      <ul>
      <li> \ref FixedPoint16Operations
      <li> \ref FixedPoint16Traits
      </ul>
  
      <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
      Namespace: vigra
  */
  template <int IntBits, FPOverflowHandling OverflowHandling>
  class FixedPoint16
  {
  public:
      static const Int32 TOTAL_BITS      = 15; // bit 16 is sign
      static const Int32 INT_BITS        = IntBits;
      static const Int32 FRACTIONAL_BITS = TOTAL_BITS - INT_BITS;
      static const Int32 MAX             = (Int32)((1u << TOTAL_BITS) - 1);
      static const Int32 MIN             = -(Int32)(1u << TOTAL_BITS);
      static const Int32 ONE             = 1 << FRACTIONAL_BITS;
      static const Int32 ONE_HALF        = ONE >> 1;
      static const Int32 FRACTIONAL_MASK = (1u << FRACTIONAL_BITS) - 1;
      static const Int32 INT_MASK        = 0xffffffffu ^ FRACTIONAL_MASK;
      
      static const FixedPoint16 zero, pi, pi_2, mpi_2; 
  
      Int16 value;
  
      FixedPoint16()
      : value(0)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Construct from an int (fractional part will become zero).
              Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
          */
      explicit FixedPoint16(Int32 v)
      : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v << FRACTIONAL_BITS))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Construct from an int by a bitwise copy. This is normally only used internally.
          */
      FixedPoint16(Int32 v, FixedPointNoShift)
      : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Construct from a double and round the fractional part to 
              <tt>FRACTIONAL_BITS</tt> accuracy. Possible overflow is handled according 
              to the target type's <tt>OverflowHandling</tt>.
          */
      explicit FixedPoint16(double rhs)
      : value(detail::FP16OverflowHandling<OverflowHandling>::exec((Int32)roundi(rhs * ONE)))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Copy constructor.
          */
      FixedPoint16(const FixedPoint16 &other)
      : value(other.value)
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Construct from a FixedPoint16 with different layout. It rounds as appropriate and 
              handles possible overflow according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2, FPOverflowHandling OverflowHandling2>
      FixedPoint16(const FixedPoint16<IntBits2, OverflowHandling2> &other)
      : value(detail::FP16OverflowHandling<OverflowHandling>::exec(
              detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value)))
      {
          VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
      }
  
          /** Assignment from int. The fractional part will become zero.  
              Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
          */
      FixedPoint16 &operator=(Int32 rhs)
      {
          value = detail::FP16OverflowHandling<OverflowHandling>::exec(rhs << FRACTIONAL_BITS);
          return *this;
      }
  
          /** Assignment form double. The fractional part is rounded, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      FixedPoint16 &operator=(double rhs)
      {
          value = detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE));
          return *this;
      }
  
          /** Copy assignment.
          */
      FixedPoint16 & operator=(const FixedPoint16 &other)
      {
          value = other.value;
          return *this;
      }
  
          /** Assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2>
      FixedPoint16 & operator=(const FixedPoint16<IntBits2, OverflowHandling> &other)
      {
          value = detail::FP16OverflowHandling<OverflowHandling>::exec(
                  detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value));
          return *this;
      }
  
          /** Conversion to float
          */
      operator float() const
      {
          return fixed_point_cast<float>(*this);
      }
  
          /** Conversion to double
          */
      operator double() const
      {
          return fixed_point_cast<double>(*this);
      }
  
          /** Unary plus.
          */
      FixedPoint16 operator+() const
      {
          return *this;
      }
  
          /** Negation.
          */
      FixedPoint16 operator-() const
      {
          return FixedPoint16(-value, FPNoShift);
      }
  
          /** Pre-increment.
          */
      FixedPoint16 & operator++()
      {
          value = detail::FP16OverflowHandling<OverflowHandling>::exec(value+ONE);
          return *this;
      }
  
          /** Post-increment.
          */
      FixedPoint16 operator++(int)
      {
          FixedPoint16 old(*this);
          ++(*this);
          return old;
      }
  
          /** Pre-decrement.
          */
      FixedPoint16 & operator--()
      {
          value = detail::FP16OverflowHandling<OverflowHandling>::exec(value-ONE);
          return *this;
      }
  
          /** Post-decrement.
          */
      FixedPoint16 operator--(int)
      {
          FixedPoint16 old(*this);
          --(*this);
          return old;
      }
  
          /** Add-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2>
      FixedPoint16 & operator+=(const FixedPoint16<IntBits2, OverflowHandling> &other)
      {
          value = detail::FP16AddImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
          return *this;
      }
  
          /** Subtract-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2>
      FixedPoint16 & operator-=(const FixedPoint16<IntBits2, OverflowHandling> &other)
      {
          value = detail::FP16SubImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
          return *this;
      }
      
          /** Multiply-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2>
      FixedPoint16 & operator*=(const FixedPoint16<IntBits2, OverflowHandling> &other)
      {
          value = detail::FP16MulImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
          return *this;
      }
      
          /** Divide-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
              handled according to the target type's <tt>OverflowHandling</tt>.
          */
      template <int IntBits2>
      FixedPoint16 & operator/=(const FixedPoint16<IntBits2, OverflowHandling> &other)
      {
          value = detail::FP16DivImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
          return *this;
      }
  };
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::zero(0);
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::pi(M_PI);
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::pi_2(0.5 * M_PI);
  
  template <int IntBits, FPOverflowHandling OverflowHandling>
  const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::mpi_2(-0.5 * M_PI); 
  
  namespace detail {
  
  template <class T>
  struct FixedPoint16Cast;
  
  #define VIGRA_FIXED_POINT_CAST(type) \
  template <> \
  struct FixedPoint16Cast<type> \
  { \
      template <int IntBits, FPOverflowHandling OverflowHandling> \
      static type cast(FixedPoint16<IntBits, OverflowHandling> v) \
      { \
          return round(v); \
      } \
  };
  
  VIGRA_FIXED_POINT_CAST(Int8)
  VIGRA_FIXED_POINT_CAST(UInt8)
  VIGRA_FIXED_POINT_CAST(Int16)
  VIGRA_FIXED_POINT_CAST(UInt16)
  VIGRA_FIXED_POINT_CAST(Int32)
  VIGRA_FIXED_POINT_CAST(UInt32)
  VIGRA_FIXED_POINT_CAST(Int64)
  VIGRA_FIXED_POINT_CAST(UInt64)
  
  #undef VIGRA_FIXED_POINT_CAST
  
  template <>
  struct FixedPoint16Cast<float>
  {
      template <int IntBits, FPOverflowHandling OverflowHandling>
      static float cast(FixedPoint16<IntBits, OverflowHandling> v)
      {
          return (float)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
      }
  };
  
  template <>
  struct FixedPoint16Cast<double>
  {
      template <int IntBits, FPOverflowHandling OverflowHandling>
      static double cast(FixedPoint16<IntBits, OverflowHandling> v)
      {
          return (double)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
      }
  };
  
  } // namespace detail
  
  /********************************************************/
  /*                                                      */
  /*                 FixedPoint16Operations               */
  /*                                                      */
  /********************************************************/
  
  /** \addtogroup FixedPoint16Operations Functions for FixedPoint16
  
      \brief     <b>\#include</b> \<vigra/fixedpoint.hxx\><br>
  
      These functions fulfill the requirements of an \ref AlgebraicRing.
  
      Namespace: vigra
      <p>
  
   */
  //@{
  
      /** Convert a FixedPoint16 to a built-in type.
          If the target is integral, the value is rounded.<br>
          Usage:
          \code
          FixedPoint16<16,15> fp(...);
          
          double d = fixed_point_cast<double>(fp);
          \endcode
      */
  template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
  TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return detail::FixedPoint16Cast<TARGET>::cast(v);
  }
  
  
      /// equal
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator==(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) == (r.value << (IntBits2 - MinIntBits));
  }
  
      /// not equal
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator!=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) != (r.value << (IntBits2 - MinIntBits));
  }
  
      /// less than
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator<(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) < (r.value << (IntBits2 - MinIntBits));
  }
  
      /// less or equal
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator<=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) <= (r.value << (IntBits2 - MinIntBits));
  }
  
      /// greater
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator>(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) > (r.value << (IntBits2 - MinIntBits));
  }
  
      /// greater or equal
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  bool operator>=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
      return (l.value << (IntBits1 - MinIntBits)) >= (r.value << (IntBits2 - MinIntBits));
  }
  
      /// addition with automatic determination of the appropriate result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
  operator+(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      typedef typename
          PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
          Result;
      return Result(detail::FP16AddImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
  }
  
      /// addition with enforced result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
  inline 
  FixedPoint16<IntBits3, OverflowHandling> &
  add(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
      FixedPoint16<IntBits3, OverflowHandling> & result)
  {
      result.value = detail::FP16AddImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
      return result;
  }
  
      /// subtraction with automatic determination of the appropriate result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
  operator-(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      typedef typename
          PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
          Result;
      return Result(detail::FP16SubImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
  }
  
      /// subtraction with enforced result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
  inline FixedPoint16<IntBits3, OverflowHandling> &
  sub(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
      FixedPoint16<IntBits3, OverflowHandling> & result)
  {
      result.value = detail::FP16SubImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
      return result;
  }
  
      /// multiplication with automatic determination of the appropriate result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
  operator*(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      typedef typename
          PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
          Result;
      return Result(detail::FP16MulImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
  }
  
      /// multiplication with enforced result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
  inline 
  FixedPoint16<IntBits3, OverflowHandling> &
  mul(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
      FixedPoint16<IntBits3, OverflowHandling> & result)
  {
      result.value = detail::FP16MulImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
      return result;
  }
  
      /// division with automatic determination of the appropriate result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
  inline
  typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
  operator/(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
  {
      typedef typename
          PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
          Result;
      return Result(detail::FP16DivImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
  }
  
      /// division with enforced result type.
  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
  inline 
  FixedPoint16<IntBits3, OverflowHandling> &
  div(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
      FixedPoint16<IntBits3, OverflowHandling> & result)
  {
      result.value = detail::FP16DivImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
      return result;
  }
  
      /// square root.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult
  sqrt(FixedPoint16<IntBits, OverflowHandling> v)
  {
      typedef typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult Result;
      enum { Shift = 15 + IntBits - 2*Result::INT_BITS };
      return Result(sqrti(v.value << Shift), FPNoShift);
  }
  
  #ifndef VIGRA_NO_HYPOT
      using ::hypot;
  #endif
  
      /// Length of hypotenuse. 
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline FixedPoint16<IntBits, OverflowHandling>
  hypot(FixedPoint16<IntBits, OverflowHandling> v1, FixedPoint16<IntBits, OverflowHandling> v2)
  {
      UInt32 l = abs(v1.value), r = abs(v2.value);   
      // sq(l) + sq(r) <= 2**31, so overflow handling after sqrti is sufficient
      return FixedPoint16<IntBits, OverflowHandling>(
                     detail::FP16OverflowHandling<OverflowHandling>::exec(sqrti(sq(l) + sq(r))), 
                     FPNoShift);
  }
  
  using std::atan2;
  
      /// Arctangent. Accuracy better than 1/3 degree (9 significant bits).
  template <int IntBits, FPOverflowHandling OverflowHandling>
  FixedPoint16<2, OverflowHandling>
  atan2(FixedPoint16<IntBits, OverflowHandling> y, FixedPoint16<IntBits, OverflowHandling> x)
  {
      enum { ResIntBits = 2 };
      typedef FixedPoint16<ResIntBits, OverflowHandling> FP;
      static const Int32 Pi_4  = 25736,       // = roundi(0.25 * M_PI * (1 << 15)), at 15 frac bits
                         Pi3_4 = 77208,       // = roundi(0.75 * M_PI * (1 << 15)),
                         c1    = 6497,        // = roundi(0.19826763260224867 * (1 << 15)), 
                         c2    = -1047730238; // = roundi(-0.9757748231899761 * (1 << 30));
                         
      // coefficients c1 and c2 minimize
      //
      // NIntegrate[(c1 r^3 + c2 r + Pi/4 - a)^4 /. r -> (Cos[a] - Sin[a])/(Cos[a] + Sin[a]), {a, 0, Pi/2}]
      //
      // Thanks to Jim Shima, http://www.dspguru.com/comp.dsp/tricks/alg/fxdatan2.htm
  
      if(x.value == 0)
          return (y.value > 0)
                     ? FP::pi_2
                     : (y.value < 0)
                           ? FP::mpi_2
                           : FP::zero;
                           
      Int32 abs_y = abs(y.value);
      Int32 r, angle;
      if(x.value > 0)
      {
          if(y.value == 0)
              return FP::zero;
          r = ((x.value - abs_y) << 15) / (x.value + abs_y); // 15 frac bits
          angle = Pi_4;
      }
      else
      {
          if(y.value == 0)
              return FP::pi;
          r = ((x.value + abs_y) << 15) / (abs_y - x.value); // 15 frac bits
          angle = Pi3_4;
      }
      
      angle += r*((c2 + c1 * (sq(r) >> 15)) >> 15) >> 15;
  
      return (y.value > 0)
                 ? FP(detail::FP16Align<0, ResIntBits, true>::exec( angle), FPNoShift)
                 : FP(detail::FP16Align<0, ResIntBits, true>::exec(-angle), FPNoShift);
  }
  
      /// absolute value.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline FixedPoint16<IntBits, OverflowHandling>
  abs(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return FixedPoint16<IntBits, OverflowHandling>(abs(v.value), FPNoShift);
  }
  
      /// squared norm (same as v*v).
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline
  typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::SquaredNormType
  squaredNorm(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return v*v;
  }
  
      /// norm (same as abs).
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline 
  typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::NormType
  norm(FixedPoint16<IntBits, OverflowHandling> const & v)
  {
      return abs(v);
  }
  
      /// fractional part. (difference between v and its floor)
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline FixedPoint16<IntBits, OverflowHandling>
  frac(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return FixedPoint16<IntBits, OverflowHandling>(
             v.value - (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK), 
             FPNoShift);
  }
  
      /// dual fractional part. (1 - frac(v))
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline FixedPoint16<IntBits, OverflowHandling>
  dual_frac(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return FixedPoint16<IntBits, OverflowHandling>(
             FixedPoint16<IntBits, OverflowHandling>::ONE - v.value + (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK), 
             FPNoShift);
  }
  
      /// rounding down.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline Int32
  floor(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return(v.value >> FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
  }
  
      /// rounding up.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline Int32
  ceil(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return((v.value + FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_MASK) >> 
                        FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
  }
  
      /// rounding to the nearest integer.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline Int32
  round(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return((v.value + FixedPoint16<IntBits, OverflowHandling>::ONE_HALF) >> 
                        FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
  }
  
      /// rounding to the nearest integer.
  template <int IntBits, FPOverflowHandling OverflowHandling>
  inline Int32
  roundi(FixedPoint16<IntBits, OverflowHandling> v)
  {
      return round(v);
  }
  
  //@}
  
  } // namespace vigra
  
  namespace std {
  
  template <int IntBits, vigra::FPOverflowHandling OverflowHandling>
  ostream & operator<<(ostream & s, vigra::FixedPoint16<IntBits, OverflowHandling> v)
  {
      s << vigra::fixed_point_cast<float>(v);
      return s;
  }
  
  } // namespace std
  
  #endif // VIGRA_FIXEDPOINT_HXX
  