[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/fixedpoint.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 2004-2005 by Ullrich Koethe                  */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036 #ifndef VIGRA_FIXEDPOINT_HXX
    00037 #define VIGRA_FIXEDPOINT_HXX
    00038
    00039 #include "mathutil.hxx"
    00040 #include "static_assert.hxx"
    00041 #include "error.hxx"
    00042 #include "numerictraits.hxx"
    00043
    00044 namespace vigra {
    00045
    00046 template <unsigned IntBits, unsigned FractionalBits>
    00047 class FixedPoint;
    00048
    00049 struct Error_FixedPointTraits_not_specialized_for_this_case;
    00050
    00051 template <class T1, class T2>
    00052 class FixedPointTraits
    00053 {
    00054 public:
    00055     typedef Error_FixedPointTraits_not_specialized_for_this_case PlusType;
    00056     typedef Error_FixedPointTraits_not_specialized_for_this_case MinusType;
    00057     typedef Error_FixedPointTraits_not_specialized_for_this_case MultipliesType;
    00058 //    typedef Error_FixedPointTraits_not_specialized_for_this_case DividesType;
    00059 };
    00060
    00061 // return type policy:
    00062 //     * try to allocate enough bits to represent the biggest possible result
    00063 //     * in case of add/subtract: if all bits of the internal int are used up,
    00064 //                                keep the representation
    00065 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00066 class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
    00067 {
    00068     enum { MaxIntBits  = (IntBits1 < IntBits2) ? IntBits2 : IntBits1,
    00069            MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1,
    00070            PlusMinusIntBits = (MaxIntBits + 1 + MaxFracBits < 32) ?
    00071                                MaxIntBits + 1 : MaxIntBits,
    00072            MultipliesFracBits = (IntBits1 + IntBits2 < 31)
    00073                                     ? (FracBits1 + FracBits2) > (31 - IntBits1 - IntBits2)
    00074                                             ? 31 - IntBits1 - IntBits2
    00075                                             : FracBits1 + FracBits2
    00076                                     : 0
    00077          };
    00078 public:
    00079     typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
    00080     typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
    00081     typedef FixedPoint<IntBits1 + IntBits2, MultipliesFracBits>  MultipliesType;
    00082 //    typedef FixedPoint<IntBits1 + FracBits2, FracBits1 + IntBits2>  DividesType;
    00083 };
    00084
    00085 template <unsigned IntBits, unsigned FracBits>
    00086 struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
    00087 {
    00088     enum { SRTotalBits = (IntBits + FracBits + 1) / 2,
    00089            SRIntBits   = (IntBits + 1) / 2,
    00090            SRFracBits  = SRTotalBits - SRIntBits
    00091          };
    00092 public:
    00093     typedef FixedPoint<IntBits, FracBits>      Type;
    00094     typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
    00095     typedef Type                               SquareRootArgument;
    00096 };
    00097
    00098
    00099 #ifndef DOXYGEN
    00100
    00101 template <int N>
    00102 struct FixedPoint_overflow_error__More_than_31_bits_requested
    00103 : staticAssert::AssertBool<(N < 32)>
    00104 {};
    00105
    00106 #endif /* DOXYGEN */
    00107
    00108
    00109
    00110 template <bool Predicate>
    00111 struct FixedPoint_assignment_error__Target_object_has_too_few_integer_bits
    00112 : staticAssert::AssertBool<Predicate>
    00113 {};
    00114
    00115 enum FixedPointNoShift { FPNoShift };
    00116
    00117 namespace detail {
    00118
    00119 template <bool MustRound>
    00120 struct FPAssignWithRound;
    00121
    00122 template <>
    00123 struct FPAssignWithRound<false>
    00124 {
    00125     template <int N>
    00126     static inline int exec(int v) { return v << (-N); }
    00127 };
    00128
    00129 template <>
    00130 struct FPAssignWithRound<true>
    00131 {
    00132     template <int N>
    00133     static inline int exec(int const v)
    00134     {
    00135         return (v + (1 << (N - 1))) >> (N);
    00136     }
    00137 };
    00138
    00139 template <bool MustRound>
    00140 struct FPMulImplementation;
    00141
    00142 template <>
    00143 struct FPMulImplementation<false>
    00144 {
    00145     template <int N>
    00146     static inline int exec(int l, int r) { return (l * r) << (-N); }
    00147 };
    00148
    00149 template <>
    00150 struct FPMulImplementation<true>
    00151 {
    00152     template <int N>
    00153     static inline int exec(int l, int r)
    00154     {
    00155         // there is not enough space in the result
    00156         // => perform calculations that preserve as much accuracy as possible
    00157         enum { diffl = N / 2, diffr = N  - diffl, maskl = (1 << diffl) - 1, maskr = (1 << diffr) - 1 };
    00158         int shiftl = l >> diffl;
    00159         int shiftr = r >> diffr;
    00160
    00161         return shiftl * shiftr + (((l & maskl) * shiftr) >> diffl) +
    00162                                  (((r & maskr) * shiftl) >> diffr);
    00163     }
    00164 };
    00165
    00166 } // namespace detail
    00167
    00168 /********************************************************/
    00169 /*                                                      */
    00170 /*                      FixedPoint                      */
    00171 /*                                                      */
    00172 /********************************************************/
    00173
    00174 /** Template for fixed point arithmetic.
    00175
    00176     Fixed point arithmetic is used when computations with fractional accuracy
    00177     must be made at the highest speed possible (e.g. in the inner loop
    00178     of a volume rendering routine). The speed-up relative to floating
    00179     point arithmetic can be dramatic, especially when one can avoid
    00180     conversions between integer and floating point numbers (these are
    00181     very expensive because integer and floating point arithmetic
    00182     resides in different pipelines).
    00183
    00184     The template wraps an <tt>int</tt> and uses <tt>IntBits</tt> to
    00185     represent the integral part of a number, and <tt>FractionalBits</tt>
    00186     for the fractional part, where <tt>IntBits + FractionalBits < 32</tt>.
    00187     (The 32rd bit is reserved because FixedPoint is a signed type).
    00188     These numbers will be automatically allocated in an intelligent way
    00189     in the result of an arithmetic operation. For example, when two
    00190     fixed point numbers are multiplied, the required number of integer
    00191     bits in the result is the sum of the number of integer bits of the
    00192     arguments, but only when so many bits are available. This is figured out
    00193     by means of FixedPointTraits, and a compile-time error is raised
    00194     when no suitable representation can be found. The idea is that the right
    00195     thing happens automatically as often as possible.
    00196
    00197     <tt>FixedPoint</tt> implements the required interface of an
    00198     \ref AlgebraicRing and the required numeric and
    00199     promotion traits. In addition, it supports functions <tt>add</tt>,
    00200     <tt>sub</tt>, and <tt>mul</tt>, where a particular layout of the result can
    00201     be enforced.
    00202
    00203     <tt>unsigned char, signed char, unsigned short, signed short, int</tt> can be
    00204     transformed into a FixedPoint with appropriate layout by means of the factory
    00205     function <tt>fixedPoint()</tt>.
    00206
    00207     <b>See also:</b>
    00208     <ul>
    00209     <li> \ref FixedPointOperations
    00210     <li> \ref FixedPointTraits
    00211     </ul>
    00212
    00213     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    00214     Namespace: vigra
    00215 */
    00216 template <unsigned IntBits, unsigned FractionalBits>
    00217 class FixedPoint
    00218 {
    00219 public:
    00220     enum {
    00221         INT_BITS        = IntBits,
    00222         FRACTIONAL_BITS = FractionalBits,
    00223         TOTAL_BITS      = IntBits + FractionalBits,
    00224         MAX             = (int)(((unsigned)1 << TOTAL_BITS) - 1),
    00225         ONE             = 1 << FractionalBits,
    00226         ONE_HALF        = ONE >> 1,
    00227         FRACTIONAL_MASK = ONE - 1,
    00228         INT_MASK        = MAX ^ FRACTIONAL_MASK
    00229     };
    00230
    00231     Int32 value;
    00232
    00233     FixedPoint()
    00234     {
    00235         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
    00236     }
    00237
    00238         /** Construct from an int (fractional part will become zero).
    00239         */
    00240     explicit FixedPoint(int v)
    00241     : value(v << FractionalBits)
    00242     {
    00243         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
    00244     }
    00245
    00246         /** Construct from an int by a bitwise copy. This is normally only used internally.
    00247         */
    00248     FixedPoint(int v, FixedPointNoShift)
    00249     : value(v)
    00250     {
    00251         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
    00252     }
    00253
    00254         /** Construct from an double and round the fractional part to
    00255             <tt>FractionalBits</tt> accuracy. A PreconditionViolation exception is raised when
    00256             the integer part is too small to represent the number.
    00257         */
    00258     explicit FixedPoint(double rhs)
    00259     : value((int)round(rhs * ONE))
    00260     {
    00261         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
    00262         vigra_precondition(abs(rhs * ONE) <= (double)MAX,
    00263             "FixedPoint(double rhs): Too few integer bits to convert rhs.");
    00264     }
    00265
    00266         /** Copy constructor.
    00267         */
    00268     FixedPoint(const FixedPoint &other)
    00269     : value(other.value)
    00270     {}
    00271
    00272         /** Construct from a FixedPoint with different layout. It rounds as appropriate and raises
    00273             a compile-time error when the target type has too few integer bits.
    00274         */
    00275     template <unsigned Int2, unsigned Frac2>
    00276     FixedPoint(const FixedPoint<Int2, Frac2> &other)
    00277     : value(detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value))
    00278     {
    00279         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
    00280         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
    00281     }
    00282
    00283         /** Assignment from int. The fractional part will become zero.
    00284             A PreconditionViolation exception is raised when
    00285             the integer part is too small to represent the number.
    00286         */
    00287     FixedPoint &operator=(int rhs)
    00288     {
    00289         vigra_precondition(abs(rhs) < (1 << IntBits),
    00290             "FixedPoint::operator=(int rhs): Too few integer bits to represent rhs.");
    00291         value = rhs << FractionalBits;
    00292         return *this;
    00293     }
    00294
    00295         /** Assignment form double. The fractional part is rounded, and a
    00296             PreconditionViolation exception is raised when
    00297             the integer part is too small to represent the number.
    00298         */
    00299     FixedPoint &operator=(double rhs)
    00300     {
    00301         vigra_precondition(abs(rhs) <= ((1 << IntBits) - 1),
    00302             "FixedPoint::operator=(double rhs): Too few integer bits to convert rhs.");
    00303         value = (int)round(rhs * ONE);
    00304         return *this;
    00305     }
    00306
    00307         /** Copy assignment.
    00308         */
    00309     FixedPoint & operator=(const FixedPoint &other)
    00310     {
    00311         value = other.value;
    00312         return *this;
    00313     }
    00314
    00315         /** Assignment from a FixedPoint with different layout. It rounds as appropriate and raises
    00316             a compile-time error when the target type has too few integer bits.
    00317         */
    00318     template <unsigned Int2, unsigned Frac2>
    00319     FixedPoint & operator=(const FixedPoint<Int2, Frac2> &other)
    00320     {
    00321         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
    00322         value = detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value);
    00323         return *this;
    00324     }
    00325
    00326         /** Negation.
    00327         */
    00328     FixedPoint operator-() const
    00329     {
    00330         return FixedPoint(-value, FPNoShift);
    00331     }
    00332
    00333         /** Pre-increment.
    00334         */
    00335     FixedPoint & operator++()
    00336     {
    00337         value += ONE;
    00338         return *this;
    00339     }
    00340
    00341         /** Post-increment.
    00342         */
    00343     FixedPoint operator++(int)
    00344     {
    00345         FixedPoint old(*this);
    00346         value += ONE;
    00347         return old;
    00348     }
    00349
    00350         /** Pre-decrement.
    00351         */
    00352     FixedPoint & operator--()
    00353     {
    00354         value -= ONE;
    00355         return *this;
    00356     }
    00357
    00358         /** Post-decrement.
    00359         */
    00360     FixedPoint operator--(int)
    00361     {
    00362         FixedPoint old(*this);
    00363         value -= ONE;
    00364         return old;
    00365     }
    00366
    00367         /** Add-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
    00368             a compile-time error when the target type has too few integer bits.
    00369         */
    00370     template <unsigned Int2, unsigned Frac2>
    00371     FixedPoint & operator+=(const FixedPoint<Int2, Frac2> &other)
    00372     {
    00373         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
    00374         value += detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
    00375         return *this;
    00376     }
    00377
    00378         /** Subtract-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
    00379             a compile-time error when the target type has too few integer bits.
    00380         */
    00381     template <unsigned Int2, unsigned Frac2>
    00382     FixedPoint & operator-=(const FixedPoint<Int2, Frac2> &other)
    00383     {
    00384         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
    00385         value -= detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
    00386         return *this;
    00387     }
    00388
    00389         /** Multiply-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
    00390             a compile-time error when the target type has too few integer bits.
    00391         */
    00392     template <unsigned Int2, unsigned Frac2>
    00393     FixedPoint & operator*=(const FixedPoint<Int2, Frac2> &other)
    00394     {
    00395         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
    00396         value = detail::FPMulImplementation<(Frac2 > 0)>::template exec<Frac2>(value, other.value);
    00397         return *this;
    00398     }
    00399 };
    00400
    00401 #define VIGRA_FIXED_POINT_FACTORY(T, INTBITS) \
    00402     inline FixedPoint<INTBITS, 0> fixedPoint(T t) \
    00403     { \
    00404         return FixedPoint<INTBITS, 0>(t, FPNoShift); \
    00405     }
    00406
    00407 VIGRA_FIXED_POINT_FACTORY(unsigned char, 8)
    00408 VIGRA_FIXED_POINT_FACTORY(signed char, 7)
    00409 VIGRA_FIXED_POINT_FACTORY(unsigned short, 16)
    00410 VIGRA_FIXED_POINT_FACTORY(signed short, 15)
    00411 VIGRA_FIXED_POINT_FACTORY(int, 31)
    00412
    00413 #undef VIGRA_FIXED_POINT_FACTORY
    00414
    00415 template <class T>
    00416 struct FixedPointCast;
    00417
    00418 #define VIGRA_FIXED_POINT_CAST(type) \
    00419 template <> \
    00420 struct FixedPointCast<type> \
    00421 { \
    00422     template <unsigned IntBits, unsigned FracBits> \
    00423     static type cast(FixedPoint<IntBits, FracBits> v) \
    00424     { \
    00425         return round(v); \
    00426     } \
    00427 };
    00428
    00429 VIGRA_FIXED_POINT_CAST(Int8)
    00430 VIGRA_FIXED_POINT_CAST(UInt8)
    00431 VIGRA_FIXED_POINT_CAST(Int16)
    00432 VIGRA_FIXED_POINT_CAST(UInt16)
    00433 VIGRA_FIXED_POINT_CAST(Int32)
    00434 VIGRA_FIXED_POINT_CAST(UInt32)
    00435
    00436 #undef VIGRA_FIXED_POINT_CAST
    00437
    00438 template <>
    00439 struct FixedPointCast<float>
    00440 {
    00441     template <unsigned IntBits, unsigned FracBits>
    00442     static float cast(FixedPoint<IntBits, FracBits> v)
    00443     {
    00444         return (float)v.value / FixedPoint<IntBits, FracBits>::ONE;
    00445     }
    00446 };
    00447
    00448 template <>
    00449 struct FixedPointCast<double>
    00450 {
    00451     template <unsigned IntBits, unsigned FracBits>
    00452     static double cast(FixedPoint<IntBits, FracBits> v)
    00453     {
    00454         return (double)v.value / FixedPoint<IntBits, FracBits>::ONE;
    00455     }
    00456 };
    00457
    00458 /********************************************************/
    00459 /*                                                      */
    00460 /*                 FixedPointOperations                 */
    00461 /*                                                      */
    00462 /********************************************************/
    00463
    00464 /** \addtogroup FixedPointOperations Functions for FixedPoint
    00465
    00466     \brief     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    00467
    00468     These functions fulfill the requirements of an \ref AlgebraicRing.
    00469
    00470     Namespace: vigra
    00471     <p>
    00472
    00473  */
    00474 //@{
    00475
    00476     /** Convert a FixedPoint to a built-in type.
    00477         If the target is integral, the value is rounded.<br>
    00478         Usage:
    00479         \code
    00480         FixedPoint<16,15> fp(...);
    00481
    00482         double d = fixed_point_cast<double>(fp);
    00483         \endcode
    00484     */
    00485 template <class TARGET, unsigned IntBits, unsigned FracBits>
    00486 TARGET fixed_point_cast(FixedPoint<IntBits, FracBits> v)
    00487 {
    00488     return FixedPointCast<TARGET>::cast(v);
    00489 }
    00490
    00491     /// equal
    00492 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00493 inline
    00494 bool operator==(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00495 {
    00496     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00497     return (l.value << (MaxFracBits - FracBits1)) == (r.value << (MaxFracBits - FracBits2));
    00498 }
    00499
    00500     /// not equal
    00501 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00502 inline
    00503 bool operator!=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00504 {
    00505     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00506     return (l.value << (MaxFracBits - FracBits1)) != (r.value << (MaxFracBits - FracBits2));
    00507 }
    00508
    00509     /// less than
    00510 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00511 inline
    00512 bool operator<(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00513 {
    00514     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00515     return (l.value << (MaxFracBits - FracBits1)) < (r.value << (MaxFracBits - FracBits2));
    00516 }
    00517
    00518     /// less or equal
    00519 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00520 inline
    00521 bool operator<=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00522 {
    00523     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00524     return (l.value << (MaxFracBits - FracBits1)) <= (r.value << (MaxFracBits - FracBits2));
    00525 }
    00526
    00527     /// greater
    00528 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00529 inline
    00530 bool operator>(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00531 {
    00532     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00533     return (l.value << (MaxFracBits - FracBits1)) > (r.value << (MaxFracBits - FracBits2));
    00534 }
    00535
    00536     /// greater or equal
    00537 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00538 inline
    00539 bool operator>=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00540 {
    00541     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00542     return (l.value << (MaxFracBits - FracBits1)) >= (r.value << (MaxFracBits - FracBits2));
    00543 }
    00544
    00545     /// addition with automatic determination of the appropriate result type.
    00546 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00547 inline
    00548 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
    00549 operator+(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00550 {
    00551     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00552     return typename
    00553         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
    00554         PlusType((l.value << (MaxFracBits - FracBits1)) + (r.value << (MaxFracBits - FracBits2)), FPNoShift);
    00555 }
    00556
    00557     /// addition with enforced result type.
    00558 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
    00559           unsigned IntBits3, unsigned FracBits3>
    00560 inline void
    00561 add(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
    00562     FixedPoint<IntBits3, FracBits3> & result)
    00563 {
    00564     result = l + r;
    00565 }
    00566
    00567     /// subtraction with automatic determination of the appropriate result type.
    00568 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00569 inline
    00570 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MinusType
    00571 operator-(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00572 {
    00573     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
    00574     return typename
    00575         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
    00576         MinusType((l.value << (MaxFracBits - FracBits1)) - (r.value << (MaxFracBits - FracBits2)), FPNoShift);
    00577 }
    00578
    00579     /// subtraction with enforced result type.
    00580 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
    00581           unsigned IntBits3, unsigned FracBits3>
    00582 inline void
    00583 sub(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
    00584     FixedPoint<IntBits3, FracBits3> & result)
    00585 {
    00586     result = l - r;
    00587 }
    00588
    00589     /// multiplication with automatic determination of the appropriate result type.
    00590 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00591 inline
    00592 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MultipliesType
    00593 operator*(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
    00594 {
    00595     typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
    00596         MultipliesType res;
    00597     mul(l, r, res);
    00598     return res;
    00599 }
    00600
    00601     /// multiplication with enforced result type.
    00602 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
    00603           unsigned IntBits3, unsigned FracBits3>
    00604 inline void
    00605 mul(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
    00606     FixedPoint<IntBits3, FracBits3> & result)
    00607 {
    00608     VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits1 + IntBits2 <= IntBits3)>));
    00609     enum { diff = FracBits1 + FracBits2 - FracBits3 };
    00610     result.value = detail::FPMulImplementation<(diff > 0)>::template exec<diff>(l.value, r.value);
    00611 }
    00612
    00613     /// square root.
    00614 template <unsigned IntBits, unsigned FracBits>
    00615 inline typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult
    00616 sqrt(FixedPoint<IntBits, FracBits> v)
    00617 {
    00618     return typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult(sqrti(v.value), FPNoShift);
    00619 }
    00620
    00621     /// absolute value.
    00622 template <unsigned IntBits, unsigned FracBits>
    00623 inline FixedPoint<IntBits, FracBits>
    00624 abs(FixedPoint<IntBits, FracBits> v)
    00625 {
    00626     return FixedPoint<IntBits, FracBits>(abs(v.value), FPNoShift);
    00627 }
    00628
    00629     /// squared norm (same as v*v).
    00630 template <unsigned IntBits, unsigned FracBits>
    00631 inline
    00632 typename FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
    00633 squaredNorm(FixedPoint<IntBits, FracBits> v)
    00634 {
    00635     return v*v;
    00636 }
    00637
    00638     /// norm (same as abs).
    00639 template <unsigned IntBits, unsigned FracBits>
    00640 inline
    00641 FixedPoint<IntBits, FracBits>
    00642 norm(FixedPoint<IntBits, FracBits> const & v)
    00643 {
    00644     return abs(v);
    00645 }
    00646
    00647     /// fractional part.
    00648 template <unsigned IntBits, unsigned FracBits>
    00649 inline FixedPoint<0, FracBits>
    00650 frac(FixedPoint<IntBits, FracBits> v)
    00651 {
    00652     return FixedPoint<0, FracBits>(v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK, FPNoShift);
    00653 }
    00654
    00655     /// dual fractional part: <tt>1 - frac(v)</tt>.
    00656 template <unsigned IntBits, unsigned FracBits>
    00657 inline FixedPoint<0, FracBits>
    00658 dual_frac(FixedPoint<IntBits, FracBits> v)
    00659 {
    00660     return FixedPoint<0, FracBits>(FixedPoint<0, FracBits>::ONE -
    00661                                    (v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK), FPNoShift);
    00662 }
    00663
    00664     /// rounding down.
    00665 template <unsigned IntBits, unsigned FracBits>
    00666 inline int
    00667 floor(FixedPoint<IntBits, FracBits> v)
    00668 {
    00669     return(v.value >> FracBits);
    00670 }
    00671
    00672     /// rounding up.
    00673 template <unsigned IntBits, unsigned FracBits>
    00674 inline int
    00675 ceil(FixedPoint<IntBits, FracBits> v)
    00676 {
    00677     return((v.value + FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK) >> FracBits);
    00678 }
    00679
    00680     /// rounding to the nearest integer.
    00681 template <unsigned IntBits, unsigned FracBits>
    00682 inline int
    00683 round(FixedPoint<IntBits, FracBits> v)
    00684 {
    00685     return((v.value + FixedPoint<IntBits, FracBits>::ONE_HALF) >> FracBits);
    00686 }
    00687
    00688 //@}
    00689
    00690 /********************************************************/
    00691 /*                                                      */
    00692 /*                     FixedPoint-Traits                */
    00693 /*                                                      */
    00694 /********************************************************/
    00695
    00696 /** \page FixedPointTraits Numeric and Promote Traits of FixedPoint
    00697
    00698     The numeric and promote traits for FixedPoint follow
    00699     the general specifications for \ref NumericPromotionTraits and
    00700     \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
    00701     partial template specialization:
    00702
    00703     \code
    00704
    00705     template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00706     class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
    00707     {
    00708         typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
    00709         typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
    00710         typedef FixedPoint<IntBits1 + IntBits2, FracBits1 + FracBits2>  MultipliesType;
    00711     };
    00712
    00713     template <unsigned IntBits, unsigned FracBits>
    00714     struct NumericTraits<FixedPoint<IntBits, FracBits> >
    00715     {
    00716         typedef FixedPoint<IntBits, FracBits> Type;
    00717             // Promote undefined because it depends on the layout, use FixedPointTraits
    00718             // RealPromote in AlgebraicRing -- multiplication with double is not supported.
    00719             // ComplexPromote in AlgebraicRing -- multiplication with double is not supported.
    00720         typedef Type ValueType;
    00721
    00722         typedef VigraFalseType isIntegral;
    00723         typedef VigraTrueType  isScalar;
    00724         typedef VigraTrueType  isSigned;
    00725         typedef VigraTrueType  isOrdered;
    00726         typedef VigraFalseType isComplex;
    00727
    00728         ... // etc.
    00729     };
    00730
    00731     template <unsigned IntBits, unsigned FracBits>
    00732     struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
    00733     {
    00734         typedef FixedPoint<IntBits, FracBits>      Type;
    00735         typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
    00736         typedef Type                               SquareRootArgument;
    00737     };
    00738
    00739     template <unsigned IntBits, unsigned FracBits>
    00740     struct NormTraits<FixedPoint<IntBits, FracBits> >
    00741     {
    00742         typedef FixedPoint<IntBits, FracBits>         Type;
    00743         typedef typename
    00744             FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
    00745                                                       SquaredNormType;
    00746         typedef Type                                  NormType;
    00747     };
    00748
    00749     template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00750     struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
    00751                          FixedPoint<IntBits2, FracBits2> >
    00752     {
    00753         typedef typename
    00754             FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
    00755             Promote;
    00756     };
    00757     \endcode
    00758
    00759     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    00760     Namespace: vigra
    00761
    00762 */
    00763 template <unsigned IntBits, unsigned FracBits>
    00764 struct NumericTraits<FixedPoint<IntBits, FracBits> >
    00765 {
    00766     typedef FixedPoint<IntBits, FracBits> Type;
    00767         //typedef FixedPoint<IntBits, FracBits> Promote;
    00768         //typedef FixedPoint<IntBits, FracBits> RealPromote;
    00769         //typedef std::complex<RealPromote> ComplexPromote;
    00770     typedef Type ValueType;
    00771
    00772     typedef VigraFalseType isIntegral;
    00773     typedef VigraTrueType  isScalar;
    00774     typedef VigraTrueType  isSigned;
    00775     typedef VigraTrueType  isOrdered;
    00776     typedef VigraFalseType isComplex;
    00777
    00778     static Type zero() { return Type(0, FPNoShift); }
    00779     static Type one() { return Type(Type::ONE, FPNoShift); }
    00780     static Type nonZero() { return one(); }
    00781     static Type epsilon() { return Type(1, FPNoShift); }
    00782     static Type smallestPositive() { return Type(1, FPNoShift); }
    00783     static Type max() { return Type( Type::MAX, FPNoShift); }
    00784     static Type min() { return -max(); }
    00785 };
    00786
    00787 template <unsigned IntBits, unsigned FracBits>
    00788 struct NormTraits<FixedPoint<IntBits, FracBits> >
    00789 {
    00790     typedef FixedPoint<IntBits, FracBits>         Type;
    00791     typedef typename
    00792         FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
    00793                                                   SquaredNormType;
    00794     typedef Type                                  NormType;
    00795 };
    00796
    00797 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
    00798 struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
    00799                      FixedPoint<IntBits2, FracBits2> >
    00800 {
    00801     typedef typename
    00802         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
    00803         Promote;
    00804 };
    00805
    00806 /***********************************************************************************/
    00807
    00808 enum FPOverflowHandling { FPOverflowIgnore, FPOverflowSaturate, FPOverflowError };
    00809
    00810 template <int IntBits, FPOverflowHandling OverflowHandling = FPOverflowIgnore>
    00811 class FixedPoint16;
    00812
    00813 /********************************************************/
    00814 /*                                                      */
    00815 /*                     FixedPoint16-Traits              */
    00816 /*                                                      */
    00817 /********************************************************/
    00818
    00819 /** \page FixedPoint16Traits Numeric and Promote Traits of FixedPoint16
    00820
    00821     The numeric and promote traits for FixedPoint16 follow
    00822     the general specifications for \ref NumericPromotionTraits and
    00823     \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
    00824     partial template specialization:
    00825
    00826     \code
    00827     template <int IntBits, FPOverflowHandling OverflowHandling>
    00828     struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
    00829     {
    00830         typedef FixedPoint16<IntBits, OverflowHandling> Type;
    00831         typedef Type                                    Promote;
    00832             // RealPromote undefined -- multiplication with double is not supported.
    00833             // ComplexPromote undefined -- multiplication with double is not supported.
    00834         typedef Type ValueType;
    00835
    00836         typedef VigraFalseType isIntegral;
    00837         typedef VigraTrueType  isScalar;
    00838         typedef VigraTrueType  isSigned;
    00839         typedef VigraTrueType  isOrdered;
    00840         typedef VigraFalseType isComplex;
    00841
    00842         ... // etc.
    00843     };
    00844
    00845     template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    00846     struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
    00847                          FixedPoint16<IntBits2, OverflowHandling> >
    00848     {
    00849         typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
    00850         ... // etc.
    00851     };
    00852
    00853     template <int IntBits, FPOverflowHandling OverflowHandling>
    00854     struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
    00855     {
    00856         typedef FixedPoint16<IntBits, OverflowHandling>     Type;
    00857         typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
    00858         typedef Type                                        NormType;
    00859     };
    00860
    00861     template <int IntBits, FPOverflowHandling OverflowHandling>
    00862     struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
    00863     {
    00864         typedef FixedPoint16<IntBits, OverflowHandling>            Type;
    00865         typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
    00866         typedef Type                                               SquareRootArgument;
    00867     };
    00868     \endcode
    00869
    00870     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    00871     Namespace: vigra
    00872
    00873 */
    00874 template <int IntBits, FPOverflowHandling OverflowHandling>
    00875 struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
    00876 {
    00877     typedef FixedPoint16<IntBits, OverflowHandling> Type;
    00878     typedef Type                                    Promote;
    00879         // RealPromote undefined -- multiplication with double is not supported.
    00880         // ComplexPromote undefined -- multiplication with double is not supported.
    00881     typedef Type ValueType;
    00882
    00883     typedef VigraFalseType isIntegral;
    00884     typedef VigraTrueType  isScalar;
    00885     typedef VigraTrueType  isSigned;
    00886     typedef VigraTrueType  isOrdered;
    00887     typedef VigraFalseType isComplex;
    00888
    00889     static Type zero() { return Type(0, FPNoShift); }
    00890     static Type one() { return Type(Type::ONE, FPNoShift); }
    00891     static Type nonZero() { return one(); }
    00892     static Type epsilon() { return Type(1, FPNoShift); }
    00893     static Type smallestPositive() { return Type(1, FPNoShift); }
    00894     static Type max() { return Type( Type::MAX, FPNoShift); }
    00895     static Type min() { return Type( Type::MIN, FPNoShift); }
    00896
    00897     static Promote toPromote(Type v) { return v; }
    00898     static Type fromPromote(Promote v) { return v; };
    00899 };
    00900
    00901 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    00902 struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
    00903                      FixedPoint16<IntBits2, OverflowHandling> >
    00904 {
    00905     typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
    00906     static Promote toPromote(FixedPoint16<IntBits1, OverflowHandling> v) { return Promote(v); }
    00907     static Promote toPromote(FixedPoint16<IntBits2, OverflowHandling> v) { return Promote(v); }
    00908 };
    00909
    00910 template <int IntBits, FPOverflowHandling OverflowHandling>
    00911 struct PromoteTraits<FixedPoint16<IntBits, OverflowHandling>,
    00912                      FixedPoint16<IntBits, OverflowHandling> >
    00913 {
    00914     typedef FixedPoint16<IntBits, OverflowHandling> Promote;
    00915     static Promote toPromote(FixedPoint16<IntBits, OverflowHandling> v) { return v; }
    00916 };
    00917
    00918 template <int IntBits, FPOverflowHandling OverflowHandling>
    00919 struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
    00920 {
    00921     typedef FixedPoint16<IntBits, OverflowHandling>     Type;
    00922     typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
    00923     typedef Type                                        NormType;
    00924 };
    00925
    00926 template <int IntBits, FPOverflowHandling OverflowHandling>
    00927 struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
    00928 {
    00929     typedef FixedPoint16<IntBits, OverflowHandling>            Type;
    00930     typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
    00931     typedef Type                                               SquareRootArgument;
    00932 };
    00933
    00934 #ifndef DOXYGEN
    00935
    00936 template <bool Compatible>
    00937 struct FixedPoint_error__Right_shift_operator_has_unsupported_semantics
    00938 : staticAssert::AssertBool<Compatible>
    00939 {};
    00940
    00941 #endif /* DOXYGEN */
    00942
    00943 template <bool Predicate>
    00944 struct FixedPoint16_assignment_error__Target_object_has_too_few_integer_bits
    00945 : staticAssert::AssertBool<Predicate>
    00946 {};
    00947
    00948 namespace detail {
    00949
    00950 template<int BeforeIntBits, int AfterIntBits,
    00951          bool Round = false,
    00952          bool RightShift = (AfterIntBits >= BeforeIntBits)>
    00953 struct FP16Align;
    00954
    00955 template<int BeforeIntBits>
    00956 struct FP16Align<BeforeIntBits, BeforeIntBits, true, true>
    00957 {
    00958     static inline Int32 exec(Int32 v)
    00959     {
    00960         return v;
    00961     }
    00962 };
    00963
    00964 template<int BeforeIntBits>
    00965 struct FP16Align<BeforeIntBits, BeforeIntBits, false, true>
    00966 {
    00967     static inline Int32 exec(Int32 v)
    00968     {
    00969         return v;
    00970     }
    00971 };
    00972
    00973 template<int BeforeIntBits, int AfterIntBits>
    00974 struct FP16Align<BeforeIntBits, AfterIntBits, false, true>
    00975 {
    00976     static inline Int32 exec(Int32 v)
    00977     {
    00978         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    00979         return v >> (AfterIntBits - BeforeIntBits);
    00980     }
    00981 };
    00982
    00983 template<int BeforeIntBits, int AfterIntBits>
    00984 struct FP16Align<BeforeIntBits, AfterIntBits, true, true>
    00985 {
    00986     enum { ONE_HALF = 1 << (AfterIntBits - BeforeIntBits - 1) };
    00987     static inline Int32 exec(Int32 v)
    00988     {
    00989         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    00990         return (v + ONE_HALF) >> (AfterIntBits - BeforeIntBits);
    00991     }
    00992 };
    00993
    00994 template<int BeforeIntBits, int AfterIntBits, bool Round>
    00995 struct FP16Align<BeforeIntBits, AfterIntBits, Round, false>
    00996 {
    00997     static inline Int32 exec(Int32 v)
    00998     {
    00999         return v << (BeforeIntBits - AfterIntBits);
    01000     }
    01001 };
    01002
    01003 template <FPOverflowHandling OverflowHandling = FPOverflowIgnore>
    01004 struct FP16OverflowHandling
    01005 {
    01006     static inline Int32 exec(Int32 v)
    01007     {
    01008         return v;
    01009     }
    01010
    01011     static inline Int32 exec(UInt32 v)
    01012     {
    01013         return v;
    01014     }
    01015 };
    01016
    01017 template <>
    01018 struct FP16OverflowHandling<FPOverflowSaturate>
    01019 {
    01020     static inline Int32 exec(Int32 v)
    01021     {
    01022         if(v >= 1 << 15)
    01023             return (1 << 15) - 1;
    01024         if(v < -(1 << 15))
    01025             return -(1 << 15);
    01026         return v;
    01027     }
    01028     static inline Int32 exec(UInt32 v)
    01029     {
    01030         if(v >= 1 << 15)
    01031             return (1 << 15) - 1;
    01032         return v;
    01033     }
    01034 };
    01035
    01036 template <>
    01037 struct FP16OverflowHandling<FPOverflowError>
    01038 {
    01039     static inline Int32 exec(Int32 v)
    01040     {
    01041         vigra_precondition(v < (1 << 15) && v >= -(1 << 15),
    01042              "FixedPoint16: Operation overflows.");
    01043         return v;
    01044     }
    01045     static inline Int32 exec(UInt32 v)
    01046     {
    01047         vigra_precondition(v < (1 << 15),
    01048              "FixedPoint16: Operation overflows.");
    01049         return v;
    01050     }
    01051 };
    01052
    01053
    01054 template <int IntBits1, int IntBits2, int IntBitsOut,
    01055           FPOverflowHandling OverflowHandling >
    01056 struct FP16AddImpl
    01057 {
    01058     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01059     static inline Int32 exec(Int32 t1, Int32 t2)
    01060     {
    01061        return FP16OverflowHandling<OverflowHandling>::exec(
    01062                   FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
    01063                       FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) +
    01064                       FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
    01065     }
    01066 };
    01067
    01068 template <int IntBits1, int IntBits2, int IntBitsOut,
    01069           FPOverflowHandling OverflowHandling >
    01070 struct FP16SubImpl
    01071 {
    01072     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01073     static inline Int32 exec(Int32 t1, Int32 t2)
    01074     {
    01075        return FP16OverflowHandling<OverflowHandling>::exec(
    01076                   FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
    01077                       FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) -
    01078                       FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
    01079     }
    01080 };
    01081
    01082 template <int IntBits1, int IntBits2, int IntBitsOut,
    01083           FPOverflowHandling OverflowHandling >
    01084 struct FP16MulImpl
    01085 {
    01086     static inline Int32 exec(Int32 t1, Int32 t2)
    01087     {
    01088         return FP16OverflowHandling<OverflowHandling>::exec(
    01089                    FP16Align<IntBits1+IntBits2, IntBitsOut+15, /*Round*/ true>::exec(t1*t2));
    01090     }
    01091 };
    01092
    01093 template <int IntBits1, int IntBits2, int IntBitsOut,
    01094           FPOverflowHandling OverflowHandling >
    01095 struct FP16DivImpl
    01096 {
    01097     static inline Int32 exec(Int32 t1, Int32 t2)
    01098     {
    01099         if(t2 == 0)
    01100             return (t1 >= 0)
    01101                        ?  (1 << 15) - 1
    01102                        : -(1 << 15);
    01103         return FP16OverflowHandling<OverflowHandling>::exec(
    01104                    FP16Align<IntBits1-IntBits2, IntBitsOut+1, /*Round*/ true>::exec((t1<<16)/t2));
    01105     }
    01106 };
    01107
    01108 } // namespace detail
    01109
    01110 /********************************************************/
    01111 /*                                                      */
    01112 /*                      FixedPoint16                    */
    01113 /*                                                      */
    01114 /********************************************************/
    01115
    01116 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
    01117 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v);
    01118
    01119 /** Template for 16-bit signed fixed point arithmetic.
    01120
    01121     Fixed point arithmetic is used when computations with fractional accuracy
    01122     must be made at the highest speed possible (e.g. in the inner loop
    01123     of a volume rendering routine). The speed-up relative to floating
    01124     point arithmetic can be dramatic, especially when one can avoid
    01125     conversions between integer and floating point numbers (these are
    01126     very expensive because integer and floating point arithmetic
    01127     resides in different pipelines).
    01128
    01129     The template wraps an <tt>Int16</tt> and uses <tt>IntBits</tt> to
    01130     represent the integral part of a number, and <tt>15 - IntBits</tt>
    01131     for the fractional part. The 16th bit is reserved because FixedPoint16
    01132     is a signed type. Results of expressions with mixed types will preserve
    01133     larger number of <tt>IntBits</tt> of the results, in order to minimize
    01134     the possibility for overflow. Nonetheless, overflow can occur, and the
    01135     template parameter <tt>OverflowHandling</tt> determines how this will be
    01136     handled:
    01137
    01138     <DL>
    01139     <DT>FPOverflowIgnore<DD> (default) Ignore overflow, i.e. use the usual modulo behavior of the
    01140                              built-in integer types.
    01141
    01142     <DT>FPOverflowSaturate<DD> Use the largest or smallest representable number (depending on sign)
    01143                                in case of overflow.
    01144
    01145     <DT>FPOverflowError<DD> Throw <tt>PreconditionViolation</tt> upon overflow. This is useful for
    01146                             debugging.
    01147     </DL>
    01148
    01149     The implementation relies on Int32-arithmetic and requires that the right-shift operator
    01150     preserves signedness. Although not enforced by the C++ standard, this is implemented
    01151     by most of today's processors. This property is checked by a
    01152     VIGRA_STATIC_ASSERT(FixedPoint_error__Right_shift_operator_has_unsupported_semantics).
    01153
    01154     <tt>FixedPoint16</tt> implements the required interface of an
    01155     \ref AlgebraicRing and the required numeric and
    01156     promotion traits. In addition, it supports functions <tt>add</tt>,
    01157     <tt>sub</tt>, <tt>mul</tt>, and <tt>div</tt>, where a particular layout
    01158     of the result can be enforced.
    01159
    01160     Built-in numeric types can be converted into <tt>FixedPoint16</tt> by the
    01161     appropriate constructors, and from <tt>FixedPoint16</tt> by means of
    01162     <tt>fixed_point_cast&lt;TargetType&gt;(fixedPoint)</tt>.
    01163
    01164     <b>See also:</b>
    01165     <ul>
    01166     <li> \ref FixedPoint16Operations
    01167     <li> \ref FixedPoint16Traits
    01168     </ul>
    01169
    01170     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    01171     Namespace: vigra
    01172 */
    01173 template <int IntBits, FPOverflowHandling OverflowHandling>
    01174 class FixedPoint16
    01175 {
    01176 public:
    01177     static const Int32 TOTAL_BITS      = 15; // bit 16 is sign
    01178     static const Int32 INT_BITS        = IntBits;
    01179     static const Int32 FRACTIONAL_BITS = TOTAL_BITS - INT_BITS;
    01180     static const Int32 MAX             = (Int32)((1u << TOTAL_BITS) - 1);
    01181     static const Int32 MIN             = -(Int32)(1u << TOTAL_BITS);
    01182     static const Int32 ONE             = 1 << FRACTIONAL_BITS;
    01183     static const Int32 ONE_HALF        = ONE >> 1;
    01184     static const Int32 FRACTIONAL_MASK = (1u << FRACTIONAL_BITS) - 1;
    01185     static const Int32 INT_MASK        = 0xffffffffu ^ FRACTIONAL_MASK;
    01186
    01187     Int16 value;
    01188
    01189     FixedPoint16()
    01190     : value(0)
    01191     {
    01192         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01193     }
    01194
    01195         /** Construct from an int (fractional part will become zero).
    01196             Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
    01197         */
    01198     explicit FixedPoint16(Int32 v)
    01199     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v << FRACTIONAL_BITS))
    01200     {
    01201         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01202     }
    01203
    01204         /** Construct from an int by a bitwise copy. This is normally only used internally.
    01205         */
    01206     FixedPoint16(Int32 v, FixedPointNoShift)
    01207     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v))
    01208     {
    01209         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01210     }
    01211
    01212         /** Construct from a double and round the fractional part to
    01213             <tt>FRACTIONAL_BITS</tt> accuracy. Possible overflow is handled according
    01214             to the target type's <tt>OverflowHandling</tt>.
    01215         */
    01216     explicit FixedPoint16(double rhs)
    01217     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE)))
    01218     {
    01219         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01220     }
    01221
    01222         /** Copy constructor.
    01223         */
    01224     FixedPoint16(const FixedPoint16 &other)
    01225     : value(other.value)
    01226     {
    01227         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01228     }
    01229
    01230         /** Construct from a FixedPoint16 with different layout. It rounds as appropriate and
    01231             handles possible overflow according to the target type's <tt>OverflowHandling</tt>.
    01232         */
    01233     template <int IntBits2, FPOverflowHandling OverflowHandling2>
    01234     FixedPoint16(const FixedPoint16<IntBits2, OverflowHandling2> &other)
    01235     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(
    01236             detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value)))
    01237     {
    01238         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
    01239     }
    01240
    01241         /** Assignment from int. The fractional part will become zero.
    01242             Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
    01243         */
    01244     FixedPoint16 &operator=(Int32 rhs)
    01245     {
    01246         value = detail::FP16OverflowHandling<OverflowHandling>::exec(rhs << FRACTIONAL_BITS);
    01247         return *this;
    01248     }
    01249
    01250         /** Assignment form double. The fractional part is rounded, and possible overflow is
    01251             handled according to the target type's <tt>OverflowHandling</tt>.
    01252         */
    01253     FixedPoint16 &operator=(double rhs)
    01254     {
    01255         value = detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE));
    01256         return *this;
    01257     }
    01258
    01259         /** Copy assignment.
    01260         */
    01261     FixedPoint16 & operator=(const FixedPoint16 &other)
    01262     {
    01263         value = other.value;
    01264         return *this;
    01265     }
    01266
    01267         /** Assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
    01268             handled according to the target type's <tt>OverflowHandling</tt>.
    01269         */
    01270     template <int IntBits2>
    01271     FixedPoint16 & operator=(const FixedPoint16<IntBits2, OverflowHandling> &other)
    01272     {
    01273         value = detail::FP16OverflowHandling<OverflowHandling>::exec(
    01274                 detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value));
    01275         return *this;
    01276     }
    01277
    01278         /** Conversion to float
    01279         */
    01280     operator float() const
    01281     {
    01282         return fixed_point_cast<float>(*this);
    01283     }
    01284
    01285         /** Conversion to double
    01286         */
    01287     operator double() const
    01288     {
    01289         return fixed_point_cast<double>(*this);
    01290     }
    01291
    01292         /** Unary plus.
    01293         */
    01294     FixedPoint16 operator+() const
    01295     {
    01296         return *this;
    01297     }
    01298
    01299         /** Negation.
    01300         */
    01301     FixedPoint16 operator-() const
    01302     {
    01303         return FixedPoint16(-value, FPNoShift);
    01304     }
    01305
    01306         /** Pre-increment.
    01307         */
    01308     FixedPoint16 & operator++()
    01309     {
    01310         value = detail::FP16OverflowHandling<OverflowHandling>::exec(value+ONE);
    01311         return *this;
    01312     }
    01313
    01314         /** Post-increment.
    01315         */
    01316     FixedPoint16 operator++(int)
    01317     {
    01318         FixedPoint16 old(*this);
    01319         ++(*this);
    01320         return old;
    01321     }
    01322
    01323         /** Pre-decrement.
    01324         */
    01325     FixedPoint16 & operator--()
    01326     {
    01327         value = detail::FP16OverflowHandling<OverflowHandling>::exec(value-ONE);
    01328         return *this;
    01329     }
    01330
    01331         /** Post-decrement.
    01332         */
    01333     FixedPoint16 operator--(int)
    01334     {
    01335         FixedPoint16 old(*this);
    01336         --(*this);
    01337         return old;
    01338     }
    01339
    01340         /** Add-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
    01341             handled according to the target type's <tt>OverflowHandling</tt>.
    01342         */
    01343     template <int IntBits2>
    01344     FixedPoint16 & operator+=(const FixedPoint16<IntBits2, OverflowHandling> &other)
    01345     {
    01346         value = detail::FP16AddImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
    01347         return *this;
    01348     }
    01349
    01350         /** Subtract-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
    01351             handled according to the target type's <tt>OverflowHandling</tt>.
    01352         */
    01353     template <int IntBits2>
    01354     FixedPoint16 & operator-=(const FixedPoint16<IntBits2, OverflowHandling> &other)
    01355     {
    01356         value = detail::FP16SubImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
    01357         return *this;
    01358     }
    01359
    01360         /** Multiply-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
    01361             handled according to the target type's <tt>OverflowHandling</tt>.
    01362         */
    01363     template <int IntBits2>
    01364     FixedPoint16 & operator*=(const FixedPoint16<IntBits2, OverflowHandling> &other)
    01365     {
    01366         value = detail::FP16MulImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
    01367         return *this;
    01368     }
    01369
    01370         /** Divide-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
    01371             handled according to the target type's <tt>OverflowHandling</tt>.
    01372         */
    01373     template <int IntBits2>
    01374     FixedPoint16 & operator/=(const FixedPoint16<IntBits2, OverflowHandling> &other)
    01375     {
    01376         value = detail::FP16DivImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
    01377         return *this;
    01378     }
    01379 };
    01380
    01381 namespace detail {
    01382
    01383 template <class T>
    01384 struct FixedPoint16Cast;
    01385
    01386 #define VIGRA_FIXED_POINT_CAST(type) \
    01387 template <> \
    01388 struct FixedPoint16Cast<type> \
    01389 { \
    01390     template <int IntBits, FPOverflowHandling OverflowHandling> \
    01391     static type cast(FixedPoint16<IntBits, OverflowHandling> v) \
    01392     { \
    01393         return round(v); \
    01394     } \
    01395 };
    01396
    01397 VIGRA_FIXED_POINT_CAST(Int8)
    01398 VIGRA_FIXED_POINT_CAST(UInt8)
    01399 VIGRA_FIXED_POINT_CAST(Int16)
    01400 VIGRA_FIXED_POINT_CAST(UInt16)
    01401 VIGRA_FIXED_POINT_CAST(Int32)
    01402 VIGRA_FIXED_POINT_CAST(UInt32)
    01403 VIGRA_FIXED_POINT_CAST(Int64)
    01404 VIGRA_FIXED_POINT_CAST(UInt64)
    01405
    01406 #undef VIGRA_FIXED_POINT_CAST
    01407
    01408 template <>
    01409 struct FixedPoint16Cast<float>
    01410 {
    01411     template <int IntBits, FPOverflowHandling OverflowHandling>
    01412     static float cast(FixedPoint16<IntBits, OverflowHandling> v)
    01413     {
    01414         return (float)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
    01415     }
    01416 };
    01417
    01418 template <>
    01419 struct FixedPoint16Cast<double>
    01420 {
    01421     template <int IntBits, FPOverflowHandling OverflowHandling>
    01422     static double cast(FixedPoint16<IntBits, OverflowHandling> v)
    01423     {
    01424         return (double)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
    01425     }
    01426 };
    01427
    01428 } // namespace detail
    01429
    01430 /********************************************************/
    01431 /*                                                      */
    01432 /*                 FixedPoint16Operations               */
    01433 /*                                                      */
    01434 /********************************************************/
    01435
    01436 /** \addtogroup FixedPoint16Operations Functions for FixedPoint16
    01437
    01438     \brief     <b>\#include</b> <vigra/fixedpoint.hxx><br>
    01439
    01440     These functions fulfill the requirements of an \ref AlgebraicRing.
    01441
    01442     Namespace: vigra
    01443     <p>
    01444
    01445  */
    01446 //@{
    01447
    01448     /** Convert a FixedPoint16 to a built-in type.
    01449         If the target is integral, the value is rounded.<br>
    01450         Usage:
    01451         \code
    01452         FixedPoint16<16,15> fp(...);
    01453
    01454         double d = fixed_point_cast<double>(fp);
    01455         \endcode
    01456     */
    01457 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
    01458 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v)
    01459 {
    01460     return detail::FixedPoint16Cast<TARGET>::cast(v);
    01461 }
    01462
    01463
    01464     /// equal
    01465 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01466 inline
    01467 bool operator==(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01468 {
    01469     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01470     return (l.value << (IntBits1 - MinIntBits)) == (r.value << (IntBits2 - MinIntBits));
    01471 }
    01472
    01473     /// not equal
    01474 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01475 inline
    01476 bool operator!=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01477 {
    01478     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01479     return (l.value << (IntBits1 - MinIntBits)) != (r.value << (IntBits2 - MinIntBits));
    01480 }
    01481
    01482     /// less than
    01483 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01484 inline
    01485 bool operator<(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01486 {
    01487     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01488     return (l.value << (IntBits1 - MinIntBits)) < (r.value << (IntBits2 - MinIntBits));
    01489 }
    01490
    01491     /// less or equal
    01492 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01493 inline
    01494 bool operator<=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01495 {
    01496     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01497     return (l.value << (IntBits1 - MinIntBits)) <= (r.value << (IntBits2 - MinIntBits));
    01498 }
    01499
    01500     /// greater
    01501 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01502 inline
    01503 bool operator>(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01504 {
    01505     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01506     return (l.value << (IntBits1 - MinIntBits)) > (r.value << (IntBits2 - MinIntBits));
    01507 }
    01508
    01509     /// greater or equal
    01510 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01511 inline
    01512 bool operator>=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01513 {
    01514     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
    01515     return (l.value << (IntBits1 - MinIntBits)) >= (r.value << (IntBits2 - MinIntBits));
    01516 }
    01517
    01518     /// addition with automatic determination of the appropriate result type.
    01519 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01520 inline
    01521 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01522 operator+(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01523 {
    01524     typedef typename
    01525         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01526         Result;
    01527     return Result(detail::FP16AddImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
    01528 }
    01529
    01530     /// addition with enforced result type.
    01531 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
    01532 inline
    01533 FixedPoint16<IntBits3, OverflowHandling> &
    01534 add(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
    01535     FixedPoint16<IntBits3, OverflowHandling> & result)
    01536 {
    01537     result.value = detail::FP16AddImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
    01538     return result;
    01539 }
    01540
    01541     /// subtraction with automatic determination of the appropriate result type.
    01542 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01543 inline
    01544 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01545 operator-(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01546 {
    01547     typedef typename
    01548         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01549         Result;
    01550     return Result(detail::FP16SubImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
    01551 }
    01552
    01553     /// subtraction with enforced result type.
    01554 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
    01555 inline FixedPoint16<IntBits3, OverflowHandling> &
    01556 sub(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
    01557     FixedPoint16<IntBits3, OverflowHandling> & result)
    01558 {
    01559     result.value = detail::FP16SubImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
    01560     return result;
    01561 }
    01562
    01563     /// multiplication with automatic determination of the appropriate result type.
    01564 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01565 inline
    01566 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01567 operator*(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01568 {
    01569     typedef typename
    01570         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01571         Result;
    01572     return Result(detail::FP16MulImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
    01573 }
    01574
    01575     /// multiplication with enforced result type.
    01576 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
    01577 inline
    01578 FixedPoint16<IntBits3, OverflowHandling> &
    01579 mul(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
    01580     FixedPoint16<IntBits3, OverflowHandling> & result)
    01581 {
    01582     result.value = detail::FP16MulImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
    01583     return result;
    01584 }
    01585
    01586     /// division with automatic determination of the appropriate result type.
    01587 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
    01588 inline
    01589 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01590 operator/(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
    01591 {
    01592     typedef typename
    01593         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
    01594         Result;
    01595     return Result(detail::FP16DivImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
    01596 }
    01597
    01598     /// division with enforced result type.
    01599 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
    01600 inline
    01601 FixedPoint16<IntBits3, OverflowHandling> &
    01602 div(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
    01603     FixedPoint16<IntBits3, OverflowHandling> & result)
    01604 {
    01605     result.value = detail::FP16DivImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
    01606     return result;
    01607 }
    01608
    01609     /// square root.
    01610 template <int IntBits, FPOverflowHandling OverflowHandling>
    01611 inline typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult
    01612 sqrt(FixedPoint16<IntBits, OverflowHandling> v)
    01613 {
    01614     typedef typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult Result;
    01615     enum { Shift = 15 + IntBits - 2*Result::INT_BITS };
    01616     return Result(sqrti(v.value << Shift), FPNoShift);
    01617 }
    01618
    01619 #ifndef VIGRA_NO_HYPOT
    01620     using ::hypot;
    01621 #endif
    01622
    01623     /// Length of hypotenuse.
    01624 template <int IntBits, FPOverflowHandling OverflowHandling>
    01625 inline FixedPoint16<IntBits, OverflowHandling>
    01626 hypot(FixedPoint16<IntBits, OverflowHandling> v1, FixedPoint16<IntBits, OverflowHandling> v2)
    01627 {
    01628     UInt32 l = abs(v1.value), r = abs(v2.value);
    01629     // sq(l) + sq(r) <= 2**31, so overflow handling after sqrti is sufficient
    01630     return FixedPoint16<IntBits, OverflowHandling>(
    01631                    detail::FP16OverflowHandling<OverflowHandling>::exec(sqrti(sq(l) + sq(r))),
    01632                    FPNoShift);
    01633 }
    01634
    01635 using std::atan2;
    01636
    01637     /// Arctangent. Accuracy better than 1/3 degree (9 significant bits).
    01638 template <int IntBits, FPOverflowHandling OverflowHandling>
    01639 FixedPoint16<2, OverflowHandling>
    01640 atan2(FixedPoint16<IntBits, OverflowHandling> y, FixedPoint16<IntBits, OverflowHandling> x)
    01641 {
    01642     enum { ResIntBits = 2 };
    01643     typedef FixedPoint16<ResIntBits, OverflowHandling> FP;
    01644     static const FP zero(0), pi(M_PI), pi_2(0.5 * M_PI), mpi_2(-0.5 * M_PI);
    01645     static const Int32 Pi_4  = roundi(0.25 * M_PI * (1 << 15)), // 15 frac bits
    01646                        Pi3_4 = roundi(0.75 * M_PI * (1 << 15)),
    01647                        c1    = roundi(0.19826763260224867 * (1 << 15)),
    01648                        c2    = roundi(-0.9757748231899761 * (1 << 30));
    01649     // coefficients c1 and c2 minimize
    01650     //
    01651     // NIntegrate[(c1 r^3 + c2 r + Pi/4 - a)^4 /. r -> (Cos[a] - Sin[a])/(Cos[a] + Sin[a]), {a, 0, Pi/2}]
    01652     //
    01653     // Thanks to Jim Shima, http://www.dspguru.com/comp.dsp/tricks/alg/fxdatan2.htm
    01654
    01655     if(x.value == 0)
    01656         return (y.value > 0)
    01657                    ? pi_2
    01658                    : (y.value < 0)
    01659                          ? mpi_2
    01660                          : zero;
    01661
    01662     Int32 abs_y = abs(y.value);
    01663     Int32 r, angle;
    01664     if(x.value > 0)
    01665     {
    01666         if(y.value == 0)
    01667             return zero;
    01668         r = ((x.value - abs_y) << 15) / (x.value + abs_y); // 15 frac bits
    01669         angle = Pi_4;
    01670     }
    01671     else
    01672     {
    01673         if(y.value == 0)
    01674             return pi;
    01675         r = ((x.value + abs_y) << 15) / (abs_y - x.value); // 15 frac bits
    01676         angle = Pi3_4;
    01677     }
    01678
    01679     angle += r*((c2 + c1 * (sq(r) >> 15)) >> 15) >> 15;
    01680
    01681     return (y.value > 0)
    01682                ? FP(detail::FP16Align<0, ResIntBits, true>::exec( angle), FPNoShift)
    01683                : FP(detail::FP16Align<0, ResIntBits, true>::exec(-angle), FPNoShift);
    01684 }
    01685
    01686     /// absolute value.
    01687 template <int IntBits, FPOverflowHandling OverflowHandling>
    01688 inline FixedPoint16<IntBits, OverflowHandling>
    01689 abs(FixedPoint16<IntBits, OverflowHandling> v)
    01690 {
    01691     return FixedPoint16<IntBits, OverflowHandling>(abs(v.value), FPNoShift);
    01692 }
    01693
    01694     /// squared norm (same as v*v).
    01695 template <int IntBits, FPOverflowHandling OverflowHandling>
    01696 inline
    01697 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::SquaredNormType
    01698 squaredNorm(FixedPoint16<IntBits, OverflowHandling> v)
    01699 {
    01700     return v*v;
    01701 }
    01702
    01703     /// norm (same as abs).
    01704 template <int IntBits, FPOverflowHandling OverflowHandling>
    01705 inline
    01706 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::NormType
    01707 norm(FixedPoint16<IntBits, OverflowHandling> const & v)
    01708 {
    01709     return abs(v);
    01710 }
    01711
    01712     /// fractional part. (difference between v and its floor)
    01713 template <int IntBits, FPOverflowHandling OverflowHandling>
    01714 inline FixedPoint16<IntBits, OverflowHandling>
    01715 frac(FixedPoint16<IntBits, OverflowHandling> v)
    01716 {
    01717     return FixedPoint16<IntBits, OverflowHandling>(
    01718            v.value - (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK),
    01719            FPNoShift);
    01720 }
    01721
    01722     /// dual fractional part. (1 - frac(v))
    01723 template <int IntBits, FPOverflowHandling OverflowHandling>
    01724 inline FixedPoint16<IntBits, OverflowHandling>
    01725 dual_frac(FixedPoint16<IntBits, OverflowHandling> v)
    01726 {
    01727     return FixedPoint16<IntBits, OverflowHandling>(
    01728            FixedPoint16<IntBits, OverflowHandling>::ONE - v.value + (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK),
    01729            FPNoShift);
    01730 }
    01731
    01732     /// rounding down.
    01733 template <int IntBits, FPOverflowHandling OverflowHandling>
    01734 inline Int32
    01735 floor(FixedPoint16<IntBits, OverflowHandling> v)
    01736 {
    01737     return(v.value >> FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
    01738 }
    01739
    01740     /// rounding up.
    01741 template <int IntBits, FPOverflowHandling OverflowHandling>
    01742 inline Int32
    01743 ceil(FixedPoint16<IntBits, OverflowHandling> v)
    01744 {
    01745     return((v.value + FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_MASK) >>
    01746                       FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
    01747 }
    01748
    01749     /// rounding to the nearest integer.
    01750 template <int IntBits, FPOverflowHandling OverflowHandling>
    01751 inline Int32
    01752 round(FixedPoint16<IntBits, OverflowHandling> v)
    01753 {
    01754     return((v.value + FixedPoint16<IntBits, OverflowHandling>::ONE_HALF) >>
    01755                       FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
    01756 }
    01757
    01758     /// rounding to the nearest integer.
    01759 template <int IntBits, FPOverflowHandling OverflowHandling>
    01760 inline Int32
    01761 roundi(FixedPoint16<IntBits, OverflowHandling> v)
    01762 {
    01763     return round(v);
    01764 }
    01765
    01766 //@}
    01767
    01768 } // namespace vigra
    01769
    01770 namespace std {
    01771
    01772 template <int IntBits, vigra::FPOverflowHandling OverflowHandling>
    01773 ostream & operator<<(ostream & s, vigra::FixedPoint16<IntBits, OverflowHandling> v)
    01774 {
    01775     s << vigra::fixed_point_cast<float>(v);
    01776     return s;
    01777 }
    01778
    01779 } // namespace std
    01780
    01781 #endif // VIGRA_FIXEDPOINT_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)