[ 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) |
html generated using doxygen and Python
|