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