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

vigra/rational.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    It was adapted from the file boost/rational.hpp of the            */
00008 /*    boost library.                                                    */
00009 /*    The VIGRA Website is                                              */
00010 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00011 /*    Please direct questions, bug reports, and contributions to        */
00012 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00013 /*        vigra@informatik.uni-hamburg.de                               */
00014 /*                                                                      */
00015 /*    Permission is hereby granted, free of charge, to any person       */
00016 /*    obtaining a copy of this software and associated documentation    */
00017 /*    files (the "Software"), to deal in the Software without           */
00018 /*    restriction, including without limitation the rights to use,      */
00019 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00020 /*    sell copies of the Software, and to permit persons to whom the    */
00021 /*    Software is furnished to do so, subject to the following          */
00022 /*    conditions:                                                       */
00023 /*                                                                      */
00024 /*    The above copyright notice and this permission notice shall be    */
00025 /*    included in all copies or substantial portions of the             */
00026 /*    Software.                                                         */
00027 /*                                                                      */
00028 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00029 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00030 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00031 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00032 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00033 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00034 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00035 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00036 /*                                                                      */
00037 /************************************************************************/
00038 
00039 // this file is based on work by Paul Moore:
00040 //
00041 //  (C) Copyright Paul Moore 1999. Permission to copy, use, modify, sell and
00042 //  distribute this software is granted provided this copyright notice appears
00043 //  in all copies. This software is provided "as is" without express or
00044 //  implied warranty, and with no claim as to its suitability for any purpose.
00045 //
00046 //  See http://www.boost.org/libs/rational for documentation.
00047 
00048 
00049 #ifndef VIGRA_RATIONAL_HPP
00050 #define VIGRA_RATIONAL_HPP
00051 
00052 #include <cmath>
00053 #include <stdexcept>
00054 #include <iosfwd>
00055 #include "config.hxx"
00056 #include "mathutil.hxx"
00057 #include "numerictraits.hxx"
00058 #include "metaprogramming.hxx"
00059 
00060 namespace vigra {
00061 
00062 /** \addtogroup MathFunctions Mathematical Functions
00063 */
00064 //@{
00065 
00066 
00067 /********************************************************/
00068 /*                                                      */
00069 /*                         gcd                          */
00070 /*                                                      */
00071 /********************************************************/
00072 
00073 /*! Calculate the greatest common divisor.
00074 
00075     This function works for arbitrary integer types, including user-defined
00076     (e.g. infinite precision) ones.
00077 
00078     <b>\#include</b> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><br>
00079     Namespace: vigra
00080 */
00081 template <typename IntType>
00082 IntType gcd(IntType n, IntType m)
00083 {
00084     // Avoid repeated construction
00085     IntType zero(0);
00086 
00087     // This is abs() - given the existence of broken compilers with Koenig
00088     // lookup issues and other problems, I code this explicitly. (Remember,
00089     // IntType may be a user-defined type).
00090     if (n < zero)
00091         n = -n;
00092     if (m < zero)
00093         m = -m;
00094 
00095     // As n and m are now positive, we can be sure that %= returns a
00096     // positive value (the standard guarantees this for built-in types,
00097     // and we require it of user-defined types).
00098     for(;;) {
00099       if(m == zero)
00100         return n;
00101       n %= m;
00102       if(n == zero)
00103         return m;
00104       m %= n;
00105     }
00106 }
00107 
00108 /********************************************************/
00109 /*                                                      */
00110 /*                         lcm                          */
00111 /*                                                      */
00112 /********************************************************/
00113 
00114 /*! Calculate the lowest common multiple.
00115 
00116     This function works for arbitrary integer types, including user-defined
00117     (e.g. infinite precision) ones.
00118 
00119     <b>\#include</b> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><br>
00120     Namespace: vigra
00121 */
00122 template <typename IntType>
00123 IntType lcm(IntType n, IntType m)
00124 {
00125     // Avoid repeated construction
00126     IntType zero(0);
00127 
00128     if (n == zero || m == zero)
00129         return zero;
00130 
00131     n /= gcd(n, m);
00132     n *= m;
00133 
00134     if (n < zero)
00135         n = -n;
00136     return n;
00137 }
00138 
00139 //@}
00140 
00141 class bad_rational : public std::domain_error
00142 {
00143 public:
00144     explicit bad_rational() : std::domain_error("bad rational: zero denominator") {}
00145 };
00146 
00147 template <typename IntType>
00148 class Rational;
00149 
00150 template <typename IntType>
00151 Rational<IntType> abs(const Rational<IntType>& r);
00152 template <typename IntType>
00153 Rational<IntType> pow(const Rational<IntType>& r, int n);
00154 template <typename IntType>
00155 Rational<IntType> floor(const Rational<IntType>& r);
00156 template <typename IntType>
00157 Rational<IntType> ceil(const Rational<IntType>& r);
00158 
00159 /********************************************************/
00160 /*                                                      */
00161 /*                       Rational                       */
00162 /*                                                      */
00163 /********************************************************/
00164 
00165 /** Template for rational numbers.
00166 
00167     This template can make use of arbitrary integer types, including
00168     user-defined (e.g. infinite precision) ones. Note, however,
00169     that overflow in either the numerator or denominator is not
00170     detected during calculations -- the standard behavior of the integer type
00171     (e.g. wrap around) applies.
00172 
00173     The class can represent and handle positive and negative infinity
00174     resulting from division by zero. Indeterminate expressions such as 0/0
00175     are signaled by a <tt>bad_rational</tt> exception which is derived from
00176     <tt>std::domain_error</tt>.
00177 
00178     <tt>Rational</tt> implements the required interface of an
00179     \ref AlgebraicField and the required \ref RationalTraits "numeric and
00180     promotion traits". All arithmetic and comparison operators, as well
00181     as the relevant algebraic functions are supported .
00182 
00183     <b>See also:</b>
00184     <ul>
00185     <li> \ref RationalTraits
00186     <li> \ref RationalOperations
00187     </ul>
00188 
00189 
00190     <b>\#include</b> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><br>
00191     Namespace: vigra
00192 */
00193 template <typename IntType>
00194 class Rational
00195 {
00196 public:
00197         /** The type of numerator and denominator
00198         */
00199     typedef IntType value_type;
00200 
00201         /** Determine whether arguments should be passed as
00202             <tt>IntType</tt> or <tt>IntType const &</tt>.
00203         */
00204     typedef typename If<typename TypeTraits<IntType>::isBuiltinType,
00205                         IntType, IntType const &>::type param_type;
00206 
00207         /** Default constructor: creates zero (<tt>0/1</tt>)
00208         */
00209     Rational()
00210     : num(0),
00211       den(1)
00212     {}
00213 
00214         /** Copy constructor
00215         */
00216     template <class U>
00217     Rational(Rational<U> const & r)
00218     : num(r.numerator()),
00219       den(r.denominator())
00220     {}
00221 
00222         /** Integer constructor: creates <tt>n/1</tt>
00223         */
00224     Rational(param_type n)
00225     : num(n),
00226       den(IntType(1))
00227     {}
00228 
00229         /** Ratio constructor: creates <tt>n/d</tt>.
00230 
00231             The ratio will be normalized unless <tt>doNormalize = false</tt>.
00232             Since the internal representation is assumed to be normalized,
00233             <tt>doNormalize = false</tt> must only be used as an optimization
00234             if <tt>n</tt> and <tt>d</tt> are known to be already normalized
00235             (i.e. have 1 as their greatest common divisor).
00236         */
00237     Rational(param_type n, param_type d, bool doNormalize = true)
00238     : num(n),
00239       den(d)
00240     {
00241         if(doNormalize)
00242             normalize();
00243     }
00244 
00245         /** Construct as an approximation of a real number.
00246 
00247             The maximal allowed relative error is given by <tt>epsilon</tt>.
00248         */
00249     explicit Rational(double v, double epsilon = 1e-4)
00250     : num(IntType(v < 0.0 ?
00251                     v/epsilon - 0.5
00252                   : v/epsilon + 0.5)),
00253       den(IntType(1.0/epsilon + 0.5))
00254     {
00255         normalize();
00256     }
00257 
00258     // Default copy constructor and assignment are fine
00259 
00260         /** Assignment from <tt>IntType</tt>.
00261         */
00262     Rational& operator=(param_type n) { return assign(n, 1); }
00263 
00264         /** Assignment from <tt>IntType</tt> pair.
00265         */
00266     Rational& assign(param_type n, param_type d, bool doNormalize = true);
00267 
00268         /** Access numerator.
00269         */
00270     param_type numerator() const { return num; }
00271 
00272         /** Access denominator.
00273         */
00274     param_type denominator() const { return den; }
00275 
00276         /** Add-assignment from <tt>Rational</tt>
00277 
00278             <tt>throws bad_rational</tt> if indeterminate expression.
00279         */
00280     Rational& operator+= (const Rational& r);
00281 
00282         /** Subtract-assignment from <tt>Rational</tt>
00283 
00284             <tt>throws bad_rational</tt> if indeterminate expression.
00285         */
00286     Rational& operator-= (const Rational& r);
00287 
00288         /** Multiply-assignment from <tt>Rational</tt>
00289 
00290             <tt>throws bad_rational</tt> if indeterminate expression.
00291         */
00292     Rational& operator*= (const Rational& r);
00293 
00294         /** Divide-assignment from <tt>Rational</tt>
00295 
00296             <tt>throws bad_rational</tt> if indeterminate expression.
00297         */
00298     Rational& operator/= (const Rational& r);
00299 
00300         /** Add-assignment from <tt>IntType</tt>
00301 
00302             <tt>throws bad_rational</tt> if indeterminate expression.
00303         */
00304     Rational& operator+= (param_type i);
00305 
00306         /** Subtract-assignment from <tt>IntType</tt>
00307 
00308             <tt>throws bad_rational</tt> if indeterminate expression.
00309         */
00310     Rational& operator-= (param_type i);
00311 
00312         /** Multiply-assignment from <tt>IntType</tt>
00313 
00314             <tt>throws bad_rational</tt> if indeterminate expression.
00315         */
00316     Rational& operator*= (param_type i);
00317 
00318         /** Divide-assignment from <tt>IntType</tt>
00319 
00320             <tt>throws bad_rational</tt> if indeterminate expression.
00321         */
00322     Rational& operator/= (param_type i);
00323 
00324         /** Pre-increment.
00325         */
00326     Rational& operator++();
00327         /** Pre-decrement.
00328         */
00329     Rational& operator--();
00330 
00331         /** Post-increment.
00332         */
00333     Rational operator++(int) { Rational res(*this); operator++(); return res; }
00334         /** Post-decrement.
00335         */
00336     Rational operator--(int) { Rational res(*this); operator--(); return res; }
00337 
00338         /** Check for zero by calling <tt>!numerator()</tt>
00339         */
00340     bool operator!() const { return !num; }
00341 
00342         /** Check whether we have positive infinity.
00343         */
00344     bool is_pinf() const
00345     {
00346         IntType zero(0);
00347         return den == zero && num > zero;
00348     }
00349 
00350         /** Check whether we have negative infinity.
00351         */
00352     bool is_ninf() const
00353     {
00354         IntType zero(0);
00355         return den == zero && num < zero;
00356     }
00357 
00358         /** Check whether we have positive or negative infinity.
00359         */
00360     bool is_inf() const
00361     {
00362         IntType zero(0);
00363         return den == zero && num != zero;
00364     }
00365 
00366         /** Check the sign.
00367 
00368             Gives 1 if the number is positive, -1 if negative, and 0 otherwise.
00369         */
00370     int sign() const
00371     {
00372         IntType zero(0);
00373         return num == zero ? 0 : num < zero ? -1 : 1;
00374     }
00375 
00376 private:
00377     // Implementation - numerator and denominator (normalized).
00378     // Other possibilities - separate whole-part, or sign, fields?
00379     IntType num;
00380     IntType den;
00381 
00382     // Representation note: Fractions are kept in normalized form at all
00383     // times. normalized form is defined as gcd(num,den) == 1 and den > 0.
00384     // In particular, note that the implementation of abs() below relies
00385     // on den always being positive.
00386     void normalize();
00387 };
00388 
00389 // Assign in place
00390 template <typename IntType>
00391 inline Rational<IntType>&
00392 Rational<IntType>::assign(param_type n, param_type d, bool doNormalize)
00393 {
00394     num = n;
00395     den = d;
00396     if(doNormalize)
00397         normalize();
00398     return *this;
00399 }
00400 
00401 // Arithmetic assignment operators
00402 template <typename IntType>
00403 Rational<IntType>& Rational<IntType>::operator+= (const Rational<IntType>& r)
00404 {
00405     IntType zero(0);
00406 
00407     // handle the Inf and NaN cases
00408     if(den == zero)
00409     {
00410         if(r.den == zero && sign()*r.sign() < 0)
00411             throw bad_rational();
00412         return *this;
00413     }
00414     if(r.den == zero)
00415     {
00416         assign(r.num, zero, false); // Inf or -Inf
00417         return *this;
00418     }
00419 
00420     // This calculation avoids overflow, and minimises the number of expensive
00421     // calculations. Thanks to Nickolay Mladenov for this algorithm.
00422     //
00423     // Proof:
00424     // We have to compute a/b + c/d, where gcd(a,b)=1 and gcd(b,c)=1.
00425     // Let g = gcd(b,d), and b = b1*g, d=d1*g. Then gcd(b1,d1)=1
00426     //
00427     // The result is (a*d1 + c*b1) / (b1*d1*g).
00428     // Now we have to normalize this ratio.
00429     // Let's assume h | gcd((a*d1 + c*b1), (b1*d1*g)), and h > 1
00430     // If h | b1 then gcd(h,d1)=1 and hence h|(a*d1+c*b1) => h|a.
00431     // But since gcd(a,b1)=1 we have h=1.
00432     // Similarly h|d1 leads to h=1.
00433     // So we have that h | gcd((a*d1 + c*b1) , (b1*d1*g)) => h|g
00434     // Finally we have gcd((a*d1 + c*b1), (b1*d1*g)) = gcd((a*d1 + c*b1), g)
00435     // Which proves that instead of normalizing the result, it is better to
00436     // divide num and den by gcd((a*d1 + c*b1), g)
00437 
00438     // Protect against self-modification
00439     IntType r_num = r.num;
00440     IntType r_den = r.den;
00441 
00442     IntType g = gcd(den, r_den);
00443     den /= g;  // = b1 from the calculations above
00444     num = num * (r_den / g) + r_num * den;
00445     g = gcd(num, g);
00446     num /= g;
00447     den *= r_den/g;
00448 
00449     return *this;
00450 }
00451 
00452 template <typename IntType>
00453 Rational<IntType>& Rational<IntType>::operator-= (const Rational<IntType>& r)
00454 {
00455     IntType zero(0);
00456 
00457     // handle the Inf and NaN cases
00458     if(den == zero)
00459     {
00460         if(r.den == zero && sign()*r.sign() > 0)
00461             throw bad_rational();
00462         return *this;
00463     }
00464     if(r.den == zero)
00465     {
00466         assign(-r.num, zero, false); // Inf or -Inf
00467         return *this;
00468     }
00469 
00470     // Protect against self-modification
00471     IntType r_num = r.num;
00472     IntType r_den = r.den;
00473 
00474     // This calculation avoids overflow, and minimises the number of expensive
00475     // calculations. It corresponds exactly to the += case above
00476     IntType g = gcd(den, r_den);
00477     den /= g;
00478     num = num * (r_den / g) - r_num * den;
00479     g = gcd(num, g);
00480     num /= g;
00481     den *= r_den/g;
00482 
00483     return *this;
00484 }
00485 
00486 template <typename IntType>
00487 Rational<IntType>& Rational<IntType>::operator*= (const Rational<IntType>& r)
00488 {
00489     IntType zero(0);
00490 
00491     // handle the Inf and NaN cases
00492     if(den == zero)
00493     {
00494         if(r.num == zero)
00495             throw bad_rational();
00496         num *= r.sign();
00497         return *this;
00498     }
00499     if(r.den == zero)
00500     {
00501         if(num == zero)
00502             throw bad_rational();
00503         num = r.num * sign();
00504         den = zero;
00505         return *this;
00506     }
00507 
00508     // Protect against self-modification
00509     IntType r_num = r.num;
00510     IntType r_den = r.den;
00511 
00512     // Avoid overflow and preserve normalization
00513     IntType gcd1 = gcd<IntType>(num, r_den);
00514     IntType gcd2 = gcd<IntType>(r_num, den);
00515     num = (num/gcd1) * (r_num/gcd2);
00516     den = (den/gcd2) * (r_den/gcd1);
00517     return *this;
00518 }
00519 
00520 template <typename IntType>
00521 Rational<IntType>& Rational<IntType>::operator/= (const Rational<IntType>& r)
00522 {
00523     IntType zero(0);
00524 
00525     // handle the Inf and NaN cases
00526     if(den == zero)
00527     {
00528         if(r.den == zero)
00529             throw bad_rational();
00530         if(r.num != zero)
00531             num *= r.sign();
00532         return *this;
00533     }
00534     if(r.num == zero)
00535     {
00536         if(num == zero)
00537             throw bad_rational();
00538         num = IntType(sign());  // normalized inf!
00539         den = zero;
00540         return *this;
00541     }
00542 
00543     if (num == zero)
00544         return *this;
00545 
00546     // Protect against self-modification
00547     IntType r_num = r.num;
00548     IntType r_den = r.den;
00549 
00550     // Avoid overflow and preserve normalization
00551     IntType gcd1 = gcd<IntType>(num, r_num);
00552     IntType gcd2 = gcd<IntType>(r_den, den);
00553     num = (num/gcd1) * (r_den/gcd2);
00554     den = (den/gcd2) * (r_num/gcd1);
00555 
00556     if (den < zero) {
00557         num = -num;
00558         den = -den;
00559     }
00560     return *this;
00561 }
00562 
00563 // Mixed-mode operators -- implement explicitly to save gcd() calculations
00564 template <typename IntType>
00565 inline Rational<IntType>&
00566 Rational<IntType>::operator+= (param_type i)
00567 {
00568     num += i * den;
00569     return *this;
00570 }
00571 
00572 template <typename IntType>
00573 inline Rational<IntType>&
00574 Rational<IntType>::operator-= (param_type i)
00575 {
00576     num -= i * den;
00577     return *this;
00578 }
00579 
00580 template <typename IntType>
00581 Rational<IntType>&
00582 Rational<IntType>::operator*= (param_type i)
00583 {
00584     if(i == IntType(1))
00585         return *this;
00586     IntType zero(0);
00587     if(i == zero)
00588     {
00589         if(den == zero)
00590         {
00591             throw bad_rational();
00592         }
00593         else
00594         {
00595             num = zero;
00596             den = IntType(1);
00597             return *this;
00598         }
00599     }
00600 
00601     IntType g = gcd(i, den);
00602     den /= g;
00603     num *= i / g;
00604     return *this;
00605 }
00606 
00607 template <typename IntType>
00608 Rational<IntType>&
00609 Rational<IntType>::operator/= (param_type i)
00610 {
00611     if(i == IntType(1))
00612         return *this;
00613 
00614     IntType zero(0);
00615     if(i == zero)
00616     {
00617         if(num == zero)
00618             throw bad_rational();
00619         num = IntType(sign()); // normalized inf!
00620         den = zero;
00621         return *this;
00622     }
00623 
00624     IntType g = gcd(i, num);
00625     if(i < zero)
00626     {
00627         num /= -g;
00628         den *= -i / g;
00629     }
00630     else
00631     {
00632         num /= g;
00633         den *= i / g;
00634     }
00635     return *this;
00636 }
00637 
00638 // Increment and decrement
00639 template <typename IntType>
00640 inline Rational<IntType>& Rational<IntType>::operator++()
00641 {
00642     // This can never denormalise the fraction
00643     num += den;
00644     return *this;
00645 }
00646 
00647 template <typename IntType>
00648 inline Rational<IntType>& Rational<IntType>::operator--()
00649 {
00650     // This can never denormalise the fraction
00651     num -= den;
00652     return *this;
00653 }
00654 
00655 // Normalisation
00656 template <typename IntType>
00657 void Rational<IntType>::normalize()
00658 {
00659     // Avoid repeated construction
00660     IntType zero(0);
00661 
00662     if (den == zero)
00663     {
00664         if(num == zero)
00665             throw bad_rational();
00666         if(num < zero)
00667             num = IntType(-1);
00668         else
00669             num = IntType(1);
00670         return;
00671     }
00672 
00673     // Handle the case of zero separately, to avoid division by zero
00674     if (num == zero) {
00675         den = IntType(1);
00676         return;
00677     }
00678 
00679     IntType g = gcd<IntType>(num, den);
00680 
00681     num /= g;
00682     den /= g;
00683 
00684     // Ensure that the denominator is positive
00685     if (den < zero) {
00686         num = -num;
00687         den = -den;
00688     }
00689 }
00690 
00691 /********************************************************/
00692 /*                                                      */
00693 /*                      Rational-Traits                 */
00694 /*                                                      */
00695 /********************************************************/
00696 
00697 /** \page RationalTraits Numeric and Promote Traits of Rational
00698 
00699     The numeric and promote traits for Rational follow
00700     the general specifications for \ref NumericPromotionTraits and
00701     \ref AlgebraicField. They are implemented in terms of the traits of the basic types by
00702     partial template specialization:
00703 
00704     \code
00705 
00706     template <class T>
00707     struct NumericTraits<Rational<T> >
00708     {
00709         typedef Rational<typename NumericTraits<T>::Promote> Promote;
00710         typedef Rational<typename NumericTraits<T>::RealPromote> RealPromote;
00711 
00712         typedef typename NumericTraits<T>::isIntegral isIntegral;
00713         typedef VigraTrueType isScalar;
00714         typedef typename NumericTraits<T>::isSigned isSigned;
00715         typedef VigraTrueType isOrdered;
00716 
00717         // etc.
00718     };
00719 
00720     template<class T>
00721     struct NormTraits<Rational<T> >
00722     {
00723         typedef Rational<T>                           Type;
00724         typedef typename NumericTraits<Type>::Promote SquaredNormType;
00725         typedef Type                                  NormType;
00726     };
00727 
00728     template <class T1, class T2>
00729     struct PromoteTraits<Rational<T1>, Rational<T2> >
00730     {
00731         typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00732     };
00733     \endcode
00734 
00735     <b>\#include</b> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><br>
00736     Namespace: vigra
00737 
00738 */
00739 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
00740 
00741 template<class T>
00742 struct NumericTraits<Rational<T> >
00743 {
00744     typedef Rational<T> Type;
00745     typedef Rational<typename NumericTraits<T>::Promote> Promote;
00746     typedef Rational<typename NumericTraits<T>::RealPromote> RealPromote;
00747     typedef std::complex<Rational<RealPromote> > ComplexPromote;
00748     typedef T ValueType;
00749 
00750     typedef typename NumericTraits<T>::isIntegral isIntegral;
00751     typedef VigraTrueType isScalar;
00752     typedef typename NumericTraits<T>::isSigned isSigned;
00753     typedef VigraTrueType isOrdered;
00754     typedef VigraFalseType isComplex;
00755 
00756     static Type zero() { return Type(0); }
00757     static Type one() { return Type(1); }
00758     static Type nonZero() { return one(); }
00759 
00760     static Promote toPromote(Type const & v)
00761         { return Promote(v.numerator(), v.denominator(), false); }
00762     static RealPromote toRealPromote(Type const & v)
00763         { return RealPromote(v.numerator(), v.denominator(), false); }
00764     static Type fromPromote(Promote const & v)
00765         { return Type(NumericTraits<T>::fromPromote(v.numerator()),
00766                       NumericTraits<T>::fromPromote(v.denominator()), false); }
00767     static Type fromRealPromote(RealPromote const & v)
00768         { return Type(NumericTraits<T>::fromRealPromote(v.numerator()),
00769                       NumericTraits<T>::fromRealPromote(v.denominator()), false); }
00770 };
00771 
00772 template<class T>
00773 struct NormTraits<Rational<T> >
00774 {
00775     typedef Rational<T>                           Type;
00776     typedef typename NumericTraits<Type>::Promote SquaredNormType;
00777     typedef Type                                  NormType;
00778 };
00779 
00780 template <class T>
00781 struct PromoteTraits<Rational<T>, Rational<T> >
00782 {
00783     typedef Rational<typename PromoteTraits<T, T>::Promote> Promote;
00784     static Promote toPromote(Rational<T> const & v) { return v; }
00785 };
00786 
00787 template <class T1, class T2>
00788 struct PromoteTraits<Rational<T1>, Rational<T2> >
00789 {
00790     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00791     static Promote toPromote(Rational<T1> const & v) { return v; }
00792     static Promote toPromote(Rational<T2> const & v) { return v; }
00793 };
00794 
00795 template <class T1, class T2>
00796 struct PromoteTraits<Rational<T1>, T2 >
00797 {
00798     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00799     static Promote toPromote(Rational<T1> const & v) { return v; }
00800     static Promote toPromote(T2 const & v) { return Promote(v); }
00801 };
00802 
00803 template <class T1, class T2>
00804 struct PromoteTraits<T1, Rational<T2> >
00805 {
00806     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00807     static Promote toPromote(T1 const & v) { return Promote(v); }
00808     static Promote toPromote(Rational<T2> const & v) { return v; }
00809 };
00810 
00811 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
00812 
00813 /********************************************************/
00814 /*                                                      */
00815 /*                   RationalOperations                 */
00816 /*                                                      */
00817 /********************************************************/
00818 
00819 /** \addtogroup RationalOperations Functions for Rational
00820 
00821     \brief     <b>\#include</b> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><br>
00822 
00823     These functions fulfill the requirements of an \ref AlgebraicField.
00824 
00825     Namespace: vigra
00826     <p>
00827 
00828  */
00829 //@{
00830 
00831 /********************************************************/
00832 /*                                                      */
00833 /*                        arithmetic                    */
00834 /*                                                      */
00835 /********************************************************/
00836 
00837     /// unary plus
00838 template <typename IntType>
00839 inline Rational<IntType> operator+ (const Rational<IntType>& r)
00840 {
00841     return r;
00842 }
00843 
00844     /// unary minus (negation)
00845 template <typename IntType>
00846 inline Rational<IntType> operator- (const Rational<IntType>& r)
00847 {
00848     return Rational<IntType>(-r.numerator(), r.denominator(), false);
00849 }
00850 
00851     /// addition
00852 template <typename IntType>
00853 inline Rational<IntType> operator+(Rational<IntType> l, Rational<IntType> const & r)
00854 {
00855     return l += r;
00856 }
00857 
00858     /// subtraction
00859 template <typename IntType>
00860 inline Rational<IntType> operator-(Rational<IntType> l, Rational<IntType> const & r)
00861 {
00862     return l -= r;
00863 }
00864 
00865     /// multiplication
00866 template <typename IntType>
00867 inline Rational<IntType> operator*(Rational<IntType> l, Rational<IntType> const & r)
00868 {
00869     return l *= r;
00870 }
00871 
00872     /// division
00873 template <typename IntType>
00874 inline Rational<IntType> operator/(Rational<IntType> l, Rational<IntType> const & r)
00875 {
00876     return l /= r;
00877 }
00878 
00879     /// addition of right-hand <tt>IntType</tt> argument
00880 template <typename IntType>
00881 inline Rational<IntType>
00882 operator+(Rational<IntType> l, typename Rational<IntType>::param_type r)
00883 {
00884     return l += r;
00885 }
00886 
00887     /// subtraction of right-hand <tt>IntType</tt> argument
00888 template <typename IntType>
00889 inline Rational<IntType>
00890 operator-(Rational<IntType> l, typename Rational<IntType>::param_type r)
00891 {
00892     return l -= r;
00893 }
00894 
00895     /// multiplication with right-hand <tt>IntType</tt> argument
00896 template <typename IntType>
00897 inline Rational<IntType>
00898 operator*(Rational<IntType> l, typename Rational<IntType>::param_type r)
00899 {
00900     return l *= r;
00901 }
00902 
00903     /// division by right-hand <tt>IntType</tt> argument
00904 template <typename IntType>
00905 inline Rational<IntType>
00906 operator/(Rational<IntType> l, typename Rational<IntType>::param_type r)
00907 {
00908     return l /= r;
00909 }
00910 
00911     /// addition of left-hand <tt>IntType</tt> argument
00912 template <typename IntType>
00913 inline Rational<IntType>
00914 operator+(typename Rational<IntType>::param_type l, Rational<IntType> r)
00915 {
00916     return r += l;
00917 }
00918 
00919     /// subtraction from left-hand <tt>IntType</tt> argument
00920 template <typename IntType>
00921 inline Rational<IntType>
00922 operator-(typename Rational<IntType>::param_type l, Rational<IntType> const & r)
00923 {
00924     return (-r) += l;
00925 }
00926 
00927     /// multiplication with left-hand <tt>IntType</tt> argument
00928 template <typename IntType>
00929 inline Rational<IntType>
00930 operator*(typename Rational<IntType>::param_type l, Rational<IntType> r)
00931 {
00932     return r *= l;
00933 }
00934 
00935     /// division of left-hand <tt>IntType</tt> argument
00936 template <typename IntType>
00937 inline Rational<IntType>
00938 operator/(typename Rational<IntType>::param_type l, Rational<IntType> const & r)
00939 {
00940     if(r.numerator() < IntType(0))
00941         return Rational<IntType>(-r.denominator(), -r.numerator(), false) *= l;
00942     else
00943         return Rational<IntType>(r.denominator(), r.numerator(), false) *= l;
00944 }
00945 
00946 /********************************************************/
00947 /*                                                      */
00948 /*                        comparison                    */
00949 /*                                                      */
00950 /********************************************************/
00951 
00952 
00953     /// equality
00954 template <typename IntType1, typename IntType2>
00955 inline bool
00956 operator== (const Rational<IntType1> & l, const Rational<IntType2>& r)
00957 {
00958     return l.denominator() == r.denominator() &&
00959            l.numerator() == r.numerator(); // works since numbers are normalized, even
00960                                            // if they represent +-infinity
00961 }
00962 
00963     /// equality with right-hand <tt>IntType2</tt> argument
00964 template <typename IntType1, typename IntType2>
00965 inline bool
00966 operator== (const Rational<IntType1> & l, IntType2 const & i)
00967 {
00968     return ((l.denominator() == IntType1(1)) && (l.numerator() == i));
00969 }
00970 
00971     /// equality with left-hand <tt>IntType1</tt> argument
00972 template <typename IntType1, typename IntType2>
00973 inline bool
00974 operator==(IntType1 const & l, Rational<IntType2> const & r)
00975 {
00976     return r == l;
00977 }
00978 
00979     /// inequality
00980 template <typename IntType1, typename IntType2>
00981 inline bool
00982 operator!=(Rational<IntType1> const & l, Rational<IntType2> const & r)
00983 {
00984     return l.denominator() != r.denominator() ||
00985            l.numerator() != r.numerator(); // works since numbers are normalized, even
00986                                            // if they represent +-infinity
00987 }
00988 
00989     /// inequality with right-hand <tt>IntType2</tt> argument
00990 template <typename IntType1, typename IntType2>
00991 inline bool
00992 operator!= (const Rational<IntType1> & l, IntType2 const & i)
00993 {
00994     return ((l.denominator() != IntType1(1)) || (l.numerator() != i));
00995 }
00996 
00997     /// inequality with left-hand <tt>IntType1</tt> argument
00998 template <typename IntType1, typename IntType2>
00999 inline bool
01000 operator!=(IntType1 const & l, Rational<IntType2> const & r)
01001 {
01002     return r != l;
01003 }
01004 
01005     /// less-than
01006 template <typename IntType1, typename IntType2>
01007 bool
01008 operator< (const Rational<IntType1> & l, const Rational<IntType2>& r)
01009 {
01010     // Avoid repeated construction
01011     typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
01012     IntType zero(0);
01013 
01014     // Handle the easy cases. Take advantage of the fact
01015     // that the denominator is never negative.
01016     if(l.denominator() == zero)
01017         if(r.denominator() == zero)
01018             // -inf < inf, !(-inf < -inf), !(inf < -inf), !(inf < inf)
01019             return l.numerator() < r.numerator();
01020         else
01021             // -inf < -1, -inf < 0, -inf < 1
01022             // !(inf < -1), !(inf < 0), !(inf < 1)
01023             return l.numerator() < zero;
01024     if(r.denominator() == zero)
01025         // -1 < inf, 0 < inf, 1 < inf
01026         // !(-1 < -inf), !(0 < -inf), !(1 < -inf)
01027         return r.numerator() > zero;
01028     // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
01029     if(l.numerator() >= zero && r.numerator() <= zero)
01030         return false;
01031     // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
01032     if(l.numerator() <= zero && r.numerator() >= zero)
01033         return true;
01034 
01035     // both numbers have the same sign (and are neither zero or +-infinity)
01036     // => calculate result, avoid overflow
01037     IntType gcd1 = gcd<IntType>(l.numerator(), r.numerator());
01038     IntType gcd2 = gcd<IntType>(r.denominator(), l.denominator());
01039     return (l.numerator()/gcd1) * (r.denominator()/gcd2) <
01040            (l.denominator()/gcd2) * (r.numerator()/gcd1);
01041 }
01042 
01043     /// less-than with right-hand <tt>IntType2</tt> argument
01044 template <typename IntType1, typename IntType2>
01045 bool
01046 operator< (const Rational<IntType1> & l, IntType2 const & i)
01047 {
01048     // Avoid repeated construction
01049     typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
01050     IntType zero(0);
01051 
01052     // Handle the easy cases. Take advantage of the fact
01053     // that the denominator is never negative.
01054     if(l.denominator() == zero)
01055         // -inf < -1, -inf < 0, -inf < 1
01056         // !(inf < -1), !(inf < 0), !(inf < 1)
01057         return l.numerator() < zero;
01058     // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
01059     if(l.numerator() >= zero && i <= zero)
01060         return false;
01061     // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
01062     if(l.numerator() <= zero && i >= zero)
01063         return true;
01064 
01065     // Now, use the fact that n/d truncates towards zero as long as n and d
01066     // are both positive.
01067     // Divide instead of multiplying to avoid overflow issues. Of course,
01068     // division may be slower, but accuracy is more important than speed...
01069     if (l.numerator() > zero)
01070         return (l.numerator()/l.denominator()) < i;
01071     else
01072         return -i < (-l.numerator()/l.denominator());
01073 }
01074 
01075     /// less-than with left-hand <tt>IntType1</tt> argument
01076 template <typename IntType1, typename IntType2>
01077 inline bool
01078 operator<(IntType1 const & l, Rational<IntType2> const & r)
01079 {
01080     return r > l;
01081 }
01082 
01083     /// greater-than
01084 template <typename IntType1, typename IntType2>
01085 inline bool
01086 operator>(Rational<IntType1> const & l, Rational<IntType2> const & r)
01087 {
01088     return r < l;
01089 }
01090 
01091     /// greater-than with right-hand <tt>IntType2</tt> argument
01092 template <typename IntType1, typename IntType2>
01093 bool
01094 operator> (const Rational<IntType1> & l, IntType2 const & i)
01095 {
01096     // Trap equality first
01097     if (l.numerator() == i && l.denominator() == IntType1(1))
01098         return false;
01099 
01100     // Otherwise, we can use operator<
01101     return !(l < i);
01102 }
01103 
01104     /// greater-than with left-hand <tt>IntType1</tt> argument
01105 template <typename IntType1, typename IntType2>
01106 inline bool
01107 operator>(IntType1 const & l, Rational<IntType2> const & r)
01108 {
01109     return r < l;
01110 }
01111 
01112     /// less-equal
01113 template <typename IntType1, typename IntType2>
01114 inline bool
01115 operator<=(Rational<IntType1> const & l, Rational<IntType2> const & r)
01116 {
01117     return !(r < l);
01118 }
01119 
01120     /// less-equal with right-hand <tt>IntType2</tt> argument
01121 template <typename IntType1, typename IntType2>
01122 inline bool
01123 operator<=(Rational<IntType1> const & l, IntType2 const & r)
01124 {
01125     return !(l > r);
01126 }
01127 
01128     /// less-equal with left-hand <tt>IntType1</tt> argument
01129 template <typename IntType1, typename IntType2>
01130 inline bool
01131 operator<=(IntType1 const & l, Rational<IntType2> const & r)
01132 {
01133     return r >= l;
01134 }
01135 
01136     /// greater-equal
01137 template <typename IntType1, typename IntType2>
01138 inline bool
01139 operator>=(Rational<IntType1> const & l, Rational<IntType2> const & r)
01140 {
01141     return !(l < r);
01142 }
01143 
01144     /// greater-equal with right-hand <tt>IntType2</tt> argument
01145 template <typename IntType1, typename IntType2>
01146 inline bool
01147 operator>=(Rational<IntType1> const & l, IntType2 const & r)
01148 {
01149     return !(l < r);
01150 }
01151 
01152     /// greater-equal with left-hand <tt>IntType1</tt> argument
01153 template <typename IntType1, typename IntType2>
01154 inline bool
01155 operator>=(IntType1 const & l, Rational<IntType2> const & r)
01156 {
01157     return r <= l;
01158 }
01159 
01160 /********************************************************/
01161 /*                                                      */
01162 /*                 algebraic functions                  */
01163 /*                                                      */
01164 /********************************************************/
01165 
01166     /// absolute value
01167 template <typename IntType>
01168 inline Rational<IntType>
01169 abs(const Rational<IntType>& r)
01170 {
01171     if (r.numerator() >= IntType(0))
01172         return r;
01173 
01174     return Rational<IntType>(-r.numerator(), r.denominator(), false);
01175 }
01176 
01177     /// norm (same as <tt>abs(r)</tt>)
01178 template <typename IntType>
01179 inline Rational<IntType>
01180 norm(const Rational<IntType>& r)
01181 {
01182     return abs(r);
01183 }
01184 
01185     /// squared norm
01186 template <typename IntType>
01187 inline typename NormTraits<Rational<IntType> >::SquaredNormType
01188 squaredNorm(const Rational<IntType>& r)
01189 {
01190     return typename NormTraits<Rational<IntType> >::SquaredNormType(sq(r.numerator()), sq(r.denominator()), false);
01191 }
01192 
01193     /** integer powers
01194 
01195         <tt>throws bad_rational</tt> if indeterminate expression.
01196     */
01197 template <typename IntType>
01198 Rational<IntType>
01199 pow(const Rational<IntType>& r, int e)
01200 {
01201     IntType zero(0);
01202     int ae;
01203     if(e == 0)
01204     {
01205         if(r.denominator() == zero)
01206                 throw bad_rational();
01207         return Rational<IntType>(IntType(1));
01208     }
01209     else if(e < 0)
01210     {
01211         if(r.numerator() == zero)
01212             return Rational<IntType>(IntType(1), zero, false);
01213         if(r.denominator() == zero)
01214             return Rational<IntType>(zero);
01215         ae = -e;
01216     }
01217     else
01218     {
01219         if(r.denominator() == zero || r.numerator() == zero)
01220             return r;
01221         ae = e;
01222     }
01223 
01224     IntType nold = r.numerator(), dold = r.denominator(),
01225             nnew = IntType(1), dnew = IntType(1);
01226     for(; ae != 0; ae >>= 1, nold *= nold, dold *= dold)
01227     {
01228         if(ae % 2 != 0)
01229         {
01230             nnew *= nold;
01231             dnew *= dold;
01232         }
01233     }
01234     if(e < 0)
01235     {
01236         if(nnew < zero)
01237             return Rational<IntType>(-dnew, -nnew, false);
01238         else
01239             return Rational<IntType>(dnew, nnew, false);
01240     }
01241     else
01242         return Rational<IntType>(nnew, dnew, false);
01243 }
01244 
01245     /// largest integer not larger than <tt>r</tt>
01246 template <typename IntType>
01247 Rational<IntType>
01248 floor(const Rational<IntType>& r)
01249 {
01250     IntType zero(0), one(1);
01251     if(r.denominator() == zero || r.denominator() == one)
01252         return r;
01253     return r.numerator() < zero ?
01254                    Rational<IntType>(r.numerator() / r.denominator() - one)
01255                  : Rational<IntType>(r.numerator() / r.denominator());
01256 }
01257 
01258     /// smallest integer not smaller than <tt>r</tt>
01259 template <typename IntType>
01260 Rational<IntType>
01261 ceil(const Rational<IntType>& r)
01262 {
01263     IntType zero(0), one(1);
01264     if(r.denominator() == zero || r.denominator() == one)
01265         return r;
01266     return r.numerator() < IntType(0) ?
01267                    Rational<IntType>(r.numerator() / r.denominator())
01268                  : Rational<IntType>(r.numerator() / r.denominator() + one);
01269 }
01270 
01271 
01272     /** Type conversion
01273 
01274         Executes <tt>static_cast<T>(numerator()) / denominator()</tt>.
01275 
01276         <b>Usage:</b>
01277 
01278         \code
01279         Rational<int> r;
01280         int i;
01281         double d;
01282         i = rational_cast<int>(r);     // round r downwards
01283         d = rational_cast<double>(r);  // represent rational as a double
01284         r = rational_cast<Rational<int> >(r);   // no change
01285         \endcode
01286     */
01287 template <typename T, typename IntType>
01288 inline T rational_cast(const Rational<IntType>& src)
01289 {
01290     return static_cast<T>(src.numerator())/src.denominator();
01291 }
01292 
01293 template <class T>
01294 inline T const & rational_cast(T const & v)
01295 {
01296     return v;
01297 }
01298 
01299 //@}
01300 
01301 template <typename IntType>
01302 std::ostream& operator<< (std::ostream& os, const vigra::Rational<IntType>& r)
01303 {
01304     os << r.numerator() << '/' << r.denominator();
01305     return os;
01306 }
01307 
01308 } // namespace vigra
01309 
01310 #endif  // VIGRA_RATIONAL_HPP
01311 

© 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.6.0 (5 Nov 2009)