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

details vigra/rational.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)