vigra - vigra/rational.hxx Source File

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

vigra/rational.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)