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

| vigra/rational.hxx |
--- | --- | ---

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)