[ 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.3, Aug 18 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 T>
00706     struct NormTraits<Rational<T> >
00707     {
00708         typedef Rational<T>                           Type;
00709         typedef typename NumericTraits<Type>::Promote SquaredNormType;
00710         typedef Type                                  NormType;
00711     };
00712 
00713     template <class T1, class T2>
00714     struct PromoteTraits<Rational<T1>, Rational<T2> >
00715     {
00716         typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00717     };
00718     \endcode
00719 
00720     <b>\#include</b> "<a href="rational_8hxx-source.html">vigra/rational.hxx</a>"<br>
00721     Namespace: vigra
00722 
00723 */
00724 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
00725 
00726 template<class T>
00727 struct NumericTraits<Rational<T> >
00728 {
00729     typedef Rational<T> Type;
00730     typedef Rational<typename NumericTraits<T>::Promote> Promote;
00731     typedef Rational<typename NumericTraits<T>::RealPromote> RealPromote;
00732     typedef std::complex<Rational<RealPromote> > ComplexPromote;
00733     typedef T ValueType;
00734 
00735     typedef typename NumericTraits<T>::isIntegral isIntegral;
00736     typedef VigraTrueType isScalar;
00737     typedef VigraTrueType isOrdered;
00738     typedef VigraFalseType isComplex;
00739 
00740     static Type zero() { return Type(0); }
00741     static Type one() { return Type(1); }
00742     static Type nonZero() { return one(); }
00743 
00744     static Promote toPromote(Type const & v)
00745         { return Promote(v.numerator(), v.denominator(), false); }
00746     static RealPromote toRealPromote(Type const & v)
00747         { return RealPromote(v.numerator(), v.denominator(), false); }
00748     static Type fromPromote(Promote const & v)
00749         { return Type(NumericTraits<T>::fromPromote(v.numerator()),
00750                       NumericTraits<T>::fromPromote(v.denominator()), false); }
00751     static Type fromRealPromote(RealPromote const & v)
00752         { return Type(NumericTraits<T>::fromRealPromote(v.numerator()),
00753                       NumericTraits<T>::fromRealPromote(v.denominator()), false); }
00754 };
00755 
00756 template<class T>
00757 struct NormTraits<Rational<T> >
00758 {
00759     typedef Rational<T>                           Type;
00760     typedef typename NumericTraits<Type>::Promote SquaredNormType;
00761     typedef Type                                  NormType;
00762 };
00763 
00764 template <class T>
00765 struct PromoteTraits<Rational<T>, Rational<T> >
00766 {
00767     typedef Rational<typename PromoteTraits<T, T>::Promote> Promote;
00768     static Promote toPromote(Rational<T> const & v) { return v; }
00769 };
00770 
00771 template <class T1, class T2>
00772 struct PromoteTraits<Rational<T1>, Rational<T2> >
00773 {
00774     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00775     static Promote toPromote(Rational<T1> const & v) { return v; }
00776     static Promote toPromote(Rational<T2> const & v) { return v; }
00777 };
00778 
00779 template <class T1, class T2>
00780 struct PromoteTraits<Rational<T1>, T2 >
00781 {
00782     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00783     static Promote toPromote(Rational<T1> const & v) { return v; }
00784     static Promote toPromote(T2 const & v) { return Promote(v); }
00785 };
00786 
00787 template <class T1, class T2>
00788 struct PromoteTraits<T1, Rational<T2> >
00789 {
00790     typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
00791     static Promote toPromote(T1 const & v) { return Promote(v); }
00792     static Promote toPromote(Rational<T2> const & v) { return v; }
00793 };
00794 
00795 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
00796 
00797 /********************************************************/
00798 /*                                                      */
00799 /*                   RationalOperations                 */
00800 /*                                                      */
00801 /********************************************************/
00802 
00803 /** \addtogroup RationalOperations Functions for Rational
00804 
00805     \brief     <b>\#include</b> "<a href="rational_8hxx-source.html">vigra/rational.hxx</a>"<br>
00806 
00807     These functions fulfill the requirements of an \ref AlgebraicField.
00808 
00809     Namespace: vigra
00810     <p>
00811 
00812  */
00813 //@{
00814 
00815 /********************************************************/
00816 /*                                                      */
00817 /*                        arithmetic                    */
00818 /*                                                      */
00819 /********************************************************/
00820 
00821     /// unary plus
00822 template <typename IntType>
00823 inline Rational<IntType> operator+ (const Rational<IntType>& r)
00824 {
00825     return r;
00826 }
00827 
00828     /// unary minus (negation)
00829 template <typename IntType>
00830 inline Rational<IntType> operator- (const Rational<IntType>& r)
00831 {
00832     return Rational<IntType>(-r.numerator(), r.denominator(), false);
00833 }
00834 
00835     /// addition
00836 template <typename IntType>
00837 inline Rational<IntType> operator+(Rational<IntType> l, Rational<IntType> const & r)
00838 {
00839     return l += r;
00840 }
00841 
00842     /// subtraction
00843 template <typename IntType>
00844 inline Rational<IntType> operator-(Rational<IntType> l, Rational<IntType> const & r)
00845 {
00846     return l -= r;
00847 }
00848 
00849     /// multiplication
00850 template <typename IntType>
00851 inline Rational<IntType> operator*(Rational<IntType> l, Rational<IntType> const & r)
00852 {
00853     return l *= r;
00854 }
00855 
00856     /// division
00857 template <typename IntType>
00858 inline Rational<IntType> operator/(Rational<IntType> l, Rational<IntType> const & r)
00859 {
00860     return l /= r;
00861 }
00862 
00863     /// addition of 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     /// subtraction of 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     /// multiplication with 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     /// division by 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     /// addition of 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     /// subtraction from 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     return (-r) += l;
00909 }
00910 
00911     /// multiplication with 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     /// division of 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     if(r.numerator() < IntType(0))
00925         return Rational<IntType>(-r.denominator(), -r.numerator(), false) *= l;
00926     else
00927         return Rational<IntType>(r.denominator(), r.numerator(), false) *= l;
00928 }
00929 
00930 /********************************************************/
00931 /*                                                      */
00932 /*                        comparison                    */
00933 /*                                                      */
00934 /********************************************************/
00935 
00936 
00937     /// equality
00938 template <typename IntType1, typename IntType2>
00939 inline bool
00940 operator== (const Rational<IntType1> & l, const Rational<IntType2>& r)
00941 {
00942     return l.denominator() == r.denominator() &&
00943            l.numerator() == r.numerator(); // works since numbers are normalized, even
00944                                            // if they represent +-infinity
00945 }
00946 
00947     /// equality with right-hand <tt>IntType2</tt> argument
00948 template <typename IntType1, typename IntType2>
00949 inline bool
00950 operator== (const Rational<IntType1> & l, IntType2 const & i)
00951 {
00952     return ((l.denominator() == IntType1(1)) && (l.numerator() == i));
00953 }
00954 
00955     /// equality with left-hand <tt>IntType1</tt> argument
00956 template <typename IntType1, typename IntType2>
00957 inline bool
00958 operator==(IntType1 const & l, Rational<IntType2> const & r)
00959 {
00960     return r == l;
00961 }
00962 
00963     /// inequality
00964 template <typename IntType1, typename IntType2>
00965 inline bool
00966 operator!=(Rational<IntType1> const & l, Rational<IntType2> const & r)
00967 {
00968     return l.denominator() != r.denominator() ||
00969            l.numerator() != r.numerator(); // works since numbers are normalized, even
00970                                            // if they represent +-infinity
00971 }
00972 
00973     /// inequality with right-hand <tt>IntType2</tt> argument
00974 template <typename IntType1, typename IntType2>
00975 inline bool
00976 operator!= (const Rational<IntType1> & l, IntType2 const & i)
00977 {
00978     return ((l.denominator() != IntType1(1)) || (l.numerator() != i));
00979 }
00980 
00981     /// inequality with left-hand <tt>IntType1</tt> argument
00982 template <typename IntType1, typename IntType2>
00983 inline bool
00984 operator!=(IntType1 const & l, Rational<IntType2> const & r)
00985 {
00986     return r != l;
00987 }
00988 
00989     /// less-than
00990 template <typename IntType1, typename IntType2>
00991 bool
00992 operator< (const Rational<IntType1> & l, const Rational<IntType2>& r)
00993 {
00994     // Avoid repeated construction
00995     typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
00996     IntType zero(0);
00997 
00998     // Handle the easy cases. Take advantage of the fact
00999     // that the denominator is never negative.
01000     if(l.denominator() == zero)
01001         if(r.denominator() == zero)
01002             // -inf < inf, !(-inf < -inf), !(inf < -inf), !(inf < inf)
01003             return l.numerator() < r.numerator();
01004         else
01005             // -inf < -1, -inf < 0, -inf < 1
01006             // !(inf < -1), !(inf < 0), !(inf < 1)
01007             return l.numerator() < zero;
01008     if(r.denominator() == zero)
01009         // -1 < inf, 0 < inf, 1 < inf
01010         // !(-1 < -inf), !(0 < -inf), !(1 < -inf)
01011         return r.numerator() > zero;
01012     // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
01013     if(l.numerator() >= zero && r.numerator() <= zero)
01014         return false;
01015     // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
01016     if(l.numerator() <= zero && r.numerator() >= zero)
01017         return true;
01018 
01019     // both numbers have the same sign (and are neither zero or +-infinity)
01020     // => calculate result, avoid overflow
01021     IntType gcd1 = gcd<IntType>(l.numerator(), r.numerator());
01022     IntType gcd2 = gcd<IntType>(r.denominator(), l.denominator());
01023     return (l.numerator()/gcd1) * (r.denominator()/gcd2) <
01024            (l.denominator()/gcd2) * (r.numerator()/gcd1);
01025 }
01026 
01027     /// less-than with right-hand <tt>IntType2</tt> argument
01028 template <typename IntType1, typename IntType2>
01029 bool
01030 operator< (const Rational<IntType1> & l, IntType2 const & i)
01031 {
01032     // Avoid repeated construction
01033     typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
01034     IntType zero(0);
01035 
01036     // Handle the easy cases. Take advantage of the fact
01037     // that the denominator is never negative.
01038     if(l.denominator() == zero)
01039         // -inf < -1, -inf < 0, -inf < 1
01040         // !(inf < -1), !(inf < 0), !(inf < 1)
01041         return l.numerator() < zero;
01042     // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
01043     if(l.numerator() >= zero && i <= zero)
01044         return false;
01045     // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
01046     if(l.numerator() <= zero && i >= zero)
01047         return true;
01048 
01049     // Now, use the fact that n/d truncates towards zero as long as n and d
01050     // are both positive.
01051     // Divide instead of multiplying to avoid overflow issues. Of course,
01052     // division may be slower, but accuracy is more important than speed...
01053     if (l.numerator() > zero)
01054         return (l.numerator()/l.denominator()) < i;
01055     else
01056         return -i < (-l.numerator()/l.denominator());
01057 }
01058 
01059     /// less-than with left-hand <tt>IntType1</tt> argument
01060 template <typename IntType1, typename IntType2>
01061 inline bool
01062 operator<(IntType1 const & l, Rational<IntType2> const & r)
01063 {
01064     return r > l;
01065 }
01066 
01067     /// greater-than
01068 template <typename IntType1, typename IntType2>
01069 inline bool
01070 operator>(Rational<IntType1> const & l, Rational<IntType2> const & r)
01071 {
01072     return r < l;
01073 }
01074 
01075     /// greater-than with right-hand <tt>IntType2</tt> argument
01076 template <typename IntType1, typename IntType2>
01077 bool
01078 operator> (const Rational<IntType1> & l, IntType2 const & i)
01079 {
01080     // Trap equality first
01081     if (l.numerator() == i && l.denominator() == IntType1(1))
01082         return false;
01083 
01084     // Otherwise, we can use operator<
01085     return !(l < i);
01086 }
01087 
01088     /// greater-than with left-hand <tt>IntType1</tt> argument
01089 template <typename IntType1, typename IntType2>
01090 inline bool
01091 operator>(IntType1 const & l, Rational<IntType2> const & r)
01092 {
01093     return r < l;
01094 }
01095 
01096     /// less-equal
01097 template <typename IntType1, typename IntType2>
01098 inline bool
01099 operator<=(Rational<IntType1> const & l, Rational<IntType2> const & r)
01100 {
01101     return !(r < l);
01102 }
01103 
01104     /// less-equal with right-hand <tt>IntType2</tt> argument
01105 template <typename IntType1, typename IntType2>
01106 inline bool
01107 operator<=(Rational<IntType1> const & l, IntType2 const & r)
01108 {
01109     return !(l > r);
01110 }
01111 
01112     /// less-equal with left-hand <tt>IntType1</tt> argument
01113 template <typename IntType1, typename IntType2>
01114 inline bool
01115 operator<=(IntType1 const & l, Rational<IntType2> const & r)
01116 {
01117     return r >= l;
01118 }
01119 
01120     /// greater-equal
01121 template <typename IntType1, typename IntType2>
01122 inline bool
01123 operator>=(Rational<IntType1> const & l, Rational<IntType2> const & r)
01124 {
01125     return !(l < r);
01126 }
01127 
01128     /// greater-equal with right-hand <tt>IntType2</tt> argument
01129 template <typename IntType1, typename IntType2>
01130 inline bool
01131 operator>=(Rational<IntType1> const & l, IntType2 const & r)
01132 {
01133     return !(l < r);
01134 }
01135 
01136     /// greater-equal with left-hand <tt>IntType1</tt> argument
01137 template <typename IntType1, typename IntType2>
01138 inline bool
01139 operator>=(IntType1 const & l, Rational<IntType2> const & r)
01140 {
01141     return r <= l;
01142 }
01143 
01144 /********************************************************/
01145 /*                                                      */
01146 /*                 algebraic functions                  */
01147 /*                                                      */
01148 /********************************************************/
01149 
01150     /// absolute value
01151 template <typename IntType>
01152 inline Rational<IntType>
01153 abs(const Rational<IntType>& r)
01154 {
01155     if (r.numerator() >= IntType(0))
01156         return r;
01157 
01158     return Rational<IntType>(-r.numerator(), r.denominator(), false);
01159 }
01160 
01161     /// norm (same as <tt>abs(r)</tt>)
01162 template <typename IntType>
01163 inline Rational<IntType>
01164 norm(const Rational<IntType>& r)
01165 {
01166     return abs(r);
01167 }
01168 
01169     /// squared norm
01170 template <typename IntType>
01171 inline typename NormTraits<Rational<IntType> >::SquaredNormType
01172 squaredNorm(const Rational<IntType>& r)
01173 {
01174     return typename NormTraits<Rational<IntType> >::SquaredNormType(sq(r.numerator()), sq(r.denominator()), false);
01175 }
01176 
01177     /** integer powers
01178 
01179         <tt>throws bad_rational</tt> if indeterminate expression.
01180     */
01181 template <typename IntType>
01182 Rational<IntType>
01183 pow(const Rational<IntType>& r, int e)
01184 {
01185     IntType zero(0);
01186     int ae;
01187     if(e == 0)
01188     {
01189         if(r.denominator() == zero)
01190                 throw bad_rational();
01191         return Rational<IntType>(IntType(1));
01192     }
01193     else if(e < 0)
01194     {
01195         if(r.numerator() == zero)
01196             return Rational<IntType>(IntType(1), zero, false);
01197         if(r.denominator() == zero)
01198             return Rational<IntType>(zero);
01199         ae = -e;
01200     }
01201     else
01202     {
01203         if(r.denominator() == zero || r.numerator() == zero)
01204             return r;
01205         ae = e;
01206     }
01207 
01208     IntType nold = r.numerator(), dold = r.denominator(),
01209             nnew = IntType(1), dnew = IntType(1);
01210     for(; ae != 0; ae >>= 1, nold *= nold, dold *= dold)
01211     {
01212         if(ae % 2 != 0)
01213         {
01214             nnew *= nold;
01215             dnew *= dold;
01216         }
01217     }
01218     if(e < 0)
01219     {
01220         if(nnew < zero)
01221             return Rational<IntType>(-dnew, -nnew, false);
01222         else
01223             return Rational<IntType>(dnew, nnew, false);
01224     }
01225     else
01226         return Rational<IntType>(nnew, dnew, false);
01227 }
01228 
01229     /// largest integer not larger than <tt>r</tt>
01230 template <typename IntType>
01231 Rational<IntType>
01232 floor(const Rational<IntType>& r)
01233 {
01234     IntType zero(0), one(1);
01235     if(r.denominator() == zero || r.denominator() == one)
01236         return r;
01237     return r.numerator() < zero ?
01238                    Rational<IntType>(r.numerator() / r.denominator() - one)
01239                  : Rational<IntType>(r.numerator() / r.denominator());
01240 }
01241 
01242     /// smallest integer not smaller than <tt>r</tt>
01243 template <typename IntType>
01244 Rational<IntType>
01245 ceil(const Rational<IntType>& r)
01246 {
01247     IntType zero(0), one(1);
01248     if(r.denominator() == zero || r.denominator() == one)
01249         return r;
01250     return r.numerator() < IntType(0) ?
01251                    Rational<IntType>(r.numerator() / r.denominator())
01252                  : Rational<IntType>(r.numerator() / r.denominator() + one);
01253 }
01254 
01255 
01256     /** Type conversion
01257 
01258         Executes <tt>static_cast<T>(numerator()) / denominator()</tt>.
01259 
01260         <b>Usage:</b>
01261 
01262         \code
01263         Rational<int> r;
01264         int i;
01265         double d;
01266         i = rational_cast<int>(r);     // round r downwards
01267         d = rational_cast<double>(r);  // represent rational as a double
01268         r = rational_cast<Rational<int> >(r);   // no change
01269         \endcode
01270     */
01271 template <typename T, typename IntType>
01272 inline T rational_cast(const Rational<IntType>& src)
01273 {
01274     return static_cast<T>(src.numerator())/src.denominator();
01275 }
01276 
01277 template <class T>
01278 inline T const & rational_cast(T const & v)
01279 {
01280     return v;
01281 }
01282 
01283 //@}
01284 
01285 } // namespace vigra
01286 
01287 template <typename IntType>
01288 std::ostream& operator<< (std::ostream& os, const vigra::Rational<IntType>& r)
01289 {
01290     os << r.numerator() << '/' << r.denominator();
01291     return os;
01292 }
01293 
01294 
01295 #endif  // VIGRA_RATIONAL_HPP
01296 

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

html generated using doxygen and Python
VIGRA 1.3.3 (18 Aug 2005)