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