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

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

© 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.7.0 (Thu Aug 25 2011)