[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><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> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><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> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><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<IntType>& 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<IntType>& 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<IntType>& 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<IntType>& 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> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><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> <<a href="rational_8hxx-source.html">vigra/rational.hxx</a>><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) |
html generated using doxygen and Python
|