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