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

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)