[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

rational.hxx
1/************************************************************************/
2/* */
3/* Copyright 1998-2004 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* It was adapted from the file boost/rational.hpp of the */
7/* boost library. */
8/* The VIGRA Website is */
9/* http://hci.iwr.uni-heidelberg.de/vigra/ */
10/* Please direct questions, bug reports, and contributions to */
11/* ullrich.koethe@iwr.uni-heidelberg.de or */
12/* vigra@informatik.uni-hamburg.de */
13/* */
14/* Permission is hereby granted, free of charge, to any person */
15/* obtaining a copy of this software and associated documentation */
16/* files (the "Software"), to deal in the Software without */
17/* restriction, including without limitation the rights to use, */
18/* copy, modify, merge, publish, distribute, sublicense, and/or */
19/* sell copies of the Software, and to permit persons to whom the */
20/* Software is furnished to do so, subject to the following */
21/* conditions: */
22/* */
23/* The above copyright notice and this permission notice shall be */
24/* included in all copies or substantial portions of the */
25/* Software. */
26/* */
27/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34/* OTHER DEALINGS IN THE SOFTWARE. */
35/* */
36/************************************************************************/
37
38// this file is based on work by Paul Moore:
39//
40// (C) Copyright Paul Moore 1999. Permission to copy, use, modify, sell and
41// distribute this software is granted provided this copyright notice appears
42// in all copies. This software is provided "as is" without express or
43// implied warranty, and with no claim as to its suitability for any purpose.
44//
45// See http://www.boost.org/libs/rational for documentation.
46
47
48#ifndef VIGRA_RATIONAL_HPP
49#define VIGRA_RATIONAL_HPP
50
51#include <cmath>
52#include <stdexcept>
53#include <iosfwd>
54#include "config.hxx"
55#include "mathutil.hxx"
56#include "numerictraits.hxx"
57#include "metaprogramming.hxx"
58
59namespace vigra {
60
61/** \addtogroup MathFunctions Mathematical Functions
62*/
63//@{
64
65
66/********************************************************/
67/* */
68/* gcd */
69/* */
70/********************************************************/
71
72/** Calculate the greatest common divisor.
73
74 This function works for arbitrary integer types, including user-defined
75 (e.g. infinite precision) ones.
76
77 <b>\#include</b> <vigra/rational.hxx><br>
78 Namespace: vigra
79*/
80template <typename IntType>
82{
83 // Avoid repeated construction
84 IntType zero(0);
85
86 // This is abs() - given the existence of broken compilers with Koenig
87 // lookup issues and other problems, I code this explicitly. (Remember,
88 // IntType may be a user-defined type).
89 if (n < zero)
90 n = -n;
91 if (m < zero)
92 m = -m;
93
94 // As n and m are now positive, we can be sure that %= returns a
95 // positive value (the standard guarantees this for built-in types,
96 // and we require it of user-defined types).
97 for(;;) {
98 if(m == zero)
99 return n;
100 n %= m;
101 if(n == zero)
102 return m;
103 m %= n;
104 }
105}
106
107/********************************************************/
108/* */
109/* lcm */
110/* */
111/********************************************************/
112
113/** Calculate the lowest common multiple.
114
115 This function works for arbitrary integer types, including user-defined
116 (e.g. infinite precision) ones.
117
118 <b>\#include</b> <vigra/rational.hxx><br>
119 Namespace: vigra
120*/
121template <typename IntType>
123{
124 // Avoid repeated construction
125 IntType zero(0);
126
127 if (n == zero || m == zero)
128 return zero;
129
130 n /= gcd(n, m);
131 n *= m;
132
133 if (n < zero)
134 n = -n;
135 return n;
136}
137
138//@}
139
140class bad_rational : public std::domain_error
141{
142public:
143 explicit bad_rational() : std::domain_error("bad rational: zero denominator") {}
144};
145
146template <typename IntType>
147class Rational;
148
149template <typename IntType>
150Rational<IntType> abs(const Rational<IntType>& r);
151template <typename IntType>
152Rational<IntType> pow(const Rational<IntType>& r, int n);
153template <typename IntType>
154Rational<IntType> floor(const Rational<IntType>& r);
155template <typename IntType>
156Rational<IntType> ceil(const Rational<IntType>& r);
157
158/********************************************************/
159/* */
160/* Rational */
161/* */
162/********************************************************/
163
164/** Template for rational numbers.
165
166 This template can make use of arbitrary integer types, including
167 user-defined (e.g. infinite precision) ones. Note, however,
168 that overflow in either the numerator or denominator is not
169 detected during calculations -- the standard behavior of the integer type
170 (e.g. wrap around) applies.
171
172 The class can represent and handle positive and negative infinity
173 resulting from division by zero. Indeterminate expressions such as 0/0
174 are signaled by a <tt>bad_rational</tt> exception which is derived from
175 <tt>std::domain_error</tt>.
176
177 <tt>Rational</tt> implements the required interface of an
178 \ref AlgebraicField and the required \ref RationalTraits "numeric and
179 promotion traits". All arithmetic and comparison operators, as well
180 as the relevant algebraic functions are supported .
181
182 <b>See also:</b>
183 <ul>
184 <li> \ref RationalTraits
185 <li> \ref RationalOperations
186 </ul>
187
188
189 <b>\#include</b> <vigra/rational.hxx><br>
190 Namespace: vigra
191*/
192template <typename IntType>
194{
195public:
196 /** The type of numerator and denominator
197 */
199
200 /** Determine whether arguments should be passed as
201 <tt>IntType</tt> or <tt>IntType const &</tt>.
202 */
204 IntType, IntType const &>::type param_type;
205
206 /** Default constructor: creates zero (<tt>0/1</tt>)
207 */
209 : num(0),
210 den(1)
211 {}
212
213 /** Copy constructor
214 */
215 template <class U>
217 : num(r.numerator()),
218 den(r.denominator())
219 {}
220
221 /** Integer constructor: creates <tt>n/1</tt>
222 */
224 : num(n),
225 den(IntType(1))
226 {}
227
228 /** Ratio constructor: creates <tt>n/d</tt>.
229
230 The ratio will be normalized unless <tt>doNormalize = false</tt>.
231 Since the internal representation is assumed to be normalized,
232 <tt>doNormalize = false</tt> must only be used as an optimization
233 if <tt>n</tt> and <tt>d</tt> are known to be already normalized
234 (i.e. have 1 as their greatest common divisor).
235 */
237 : num(n),
238 den(d)
239 {
240 if(doNormalize)
241 normalize();
242 }
243
244 /** Construct as an approximation of a real number.
245
246 The maximal allowed relative error is given by <tt>epsilon</tt>.
247 */
248 explicit Rational(double v, double epsilon = 1e-4)
249 : num(IntType(v < 0.0 ?
250 v/epsilon - 0.5
251 : v/epsilon + 0.5)),
252 den(IntType(1.0/epsilon + 0.5))
253 {
254 normalize();
255 }
256
257 // Default copy constructor and assignment are fine
258
259 /** Assignment from <tt>IntType</tt>.
260 */
261 Rational& operator=(param_type n) { return assign(n, 1); }
262
263 /** Assignment from <tt>IntType</tt> pair.
264 */
265 Rational& assign(param_type n, param_type d, bool doNormalize = true);
266
267 /** Access numerator.
268 */
269 param_type numerator() const { return num; }
270
271 /** Access denominator.
272 */
273 param_type denominator() const { return den; }
274
275 /** Add-assignment from <tt>Rational</tt>
276
277 <tt>throws bad_rational</tt> if indeterminate expression.
278 */
279 Rational& operator+= (const Rational & r);
280
281 /** Subtract-assignment from <tt>Rational</tt>
282
283 <tt>throws bad_rational</tt> if indeterminate expression.
284 */
285 Rational& operator-= (const Rational & r);
286
287 /** Multiply-assignment from <tt>Rational</tt>
288
289 <tt>throws bad_rational</tt> if indeterminate expression.
290 */
291 Rational& operator*= (const Rational & r);
292
293 /** Divide-assignment from <tt>Rational</tt>
294
295 <tt>throws bad_rational</tt> if indeterminate expression.
296 */
297 Rational& operator/= (const Rational & r);
298
299 /** Add-assignment from <tt>IntType</tt>
300
301 <tt>throws bad_rational</tt> if indeterminate expression.
302 */
304
305 /** Subtract-assignment from <tt>IntType</tt>
306
307 <tt>throws bad_rational</tt> if indeterminate expression.
308 */
310
311 /** Multiply-assignment from <tt>IntType</tt>
312
313 <tt>throws bad_rational</tt> if indeterminate expression.
314 */
316
317 /** Divide-assignment from <tt>IntType</tt>
318
319 <tt>throws bad_rational</tt> if indeterminate expression.
320 */
322
323 /** Pre-increment.
324 */
326 /** Pre-decrement.
327 */
329
330 /** Post-increment.
331 */
332 Rational operator++(int) { Rational res(*this); operator++(); return res; }
333 /** Post-decrement.
334 */
335 Rational operator--(int) { Rational res(*this); operator--(); return res; }
336
337 /** Check for zero by calling <tt>!numerator()</tt>
338 */
339 bool operator!() const { return !num; }
340
341 /** Check whether we have positive infinity.
342 */
343 bool is_pinf() const
344 {
345 IntType zero(0);
346 return den == zero && num > zero;
347 }
348
349 /** Check whether we have negative infinity.
350 */
351 bool is_ninf() const
352 {
353 IntType zero(0);
354 return den == zero && num < zero;
355 }
356
357 /** Check whether we have positive or negative infinity.
358 */
359 bool is_inf() const
360 {
361 IntType zero(0);
362 return den == zero && num != zero;
363 }
364
365 /** Check the sign.
366
367 Gives 1 if the number is positive, -1 if negative, and 0 otherwise.
368 */
369 int sign() const
370 {
371 IntType zero(0);
372 return num == zero ? 0 : num < zero ? -1 : 1;
373 }
374
375private:
376 // Implementation - numerator and denominator (normalized).
377 // Other possibilities - separate whole-part, or sign, fields?
378 IntType num;
379 IntType den;
380
381 // Representation note: Fractions are kept in normalized form at all
382 // times. normalized form is defined as gcd(num,den) == 1 and den > 0.
383 // In particular, note that the implementation of abs() below relies
384 // on den always being positive.
385 void normalize();
386};
387
388// Assign in place
389template <typename IntType>
390inline Rational<IntType>&
392{
393 num = n;
394 den = d;
395 if(doNormalize)
396 normalize();
397 return *this;
398}
399
400// Arithmetic assignment operators
401template <typename IntType>
403{
404 IntType zero(0);
405
406 // handle the Inf and NaN cases
407 if(den == zero)
408 {
409 if(r.den == zero && sign()*r.sign() < 0)
410 throw bad_rational();
411 return *this;
412 }
413 if(r.den == zero)
414 {
415 assign(r.num, zero, false); // Inf or -Inf
416 return *this;
417 }
418
419 // This calculation avoids overflow, and minimises the number of expensive
420 // calculations. Thanks to Nickolay Mladenov for this algorithm.
421 //
422 // Proof:
423 // We have to compute a/b + c/d, where gcd(a,b)=1 and gcd(b,c)=1.
424 // Let g = gcd(b,d), and b = b1*g, d=d1*g. Then gcd(b1,d1)=1
425 //
426 // The result is (a*d1 + c*b1) / (b1*d1*g).
427 // Now we have to normalize this ratio.
428 // Let's assume h | gcd((a*d1 + c*b1), (b1*d1*g)), and h > 1
429 // If h | b1 then gcd(h,d1)=1 and hence h|(a*d1+c*b1) => h|a.
430 // But since gcd(a,b1)=1 we have h=1.
431 // Similarly h|d1 leads to h=1.
432 // So we have that h | gcd((a*d1 + c*b1) , (b1*d1*g)) => h|g
433 // Finally we have gcd((a*d1 + c*b1), (b1*d1*g)) = gcd((a*d1 + c*b1), g)
434 // Which proves that instead of normalizing the result, it is better to
435 // divide num and den by gcd((a*d1 + c*b1), g)
436
437 // Protect against self-modification
438 IntType r_num = r.num;
439 IntType r_den = r.den;
440
441 IntType g = gcd(den, r_den);
442 den /= g; // = b1 from the calculations above
443 num = num * (r_den / g) + r_num * den;
444 g = gcd(num, g);
445 num /= g;
446 den *= r_den/g;
447
448 return *this;
449}
450
451template <typename IntType>
453{
454 IntType zero(0);
455
456 // handle the Inf and NaN cases
457 if(den == zero)
458 {
459 if(r.den == zero && sign()*r.sign() > 0)
460 throw bad_rational();
461 return *this;
462 }
463 if(r.den == zero)
464 {
465 assign(-r.num, zero, false); // Inf or -Inf
466 return *this;
467 }
468
469 // Protect against self-modification
470 IntType r_num = r.num;
471 IntType r_den = r.den;
472
473 // This calculation avoids overflow, and minimises the number of expensive
474 // calculations. It corresponds exactly to the += case above
475 IntType g = gcd(den, r_den);
476 den /= g;
477 num = num * (r_den / g) - r_num * den;
478 g = gcd(num, g);
479 num /= g;
480 den *= r_den/g;
481
482 return *this;
483}
484
485template <typename IntType>
487{
488 IntType zero(0);
489
490 // handle the Inf and NaN cases
491 if(den == zero)
492 {
493 if(r.num == zero)
494 throw bad_rational();
495 num *= r.sign();
496 return *this;
497 }
498 if(r.den == zero)
499 {
500 if(num == zero)
501 throw bad_rational();
502 num = r.num * sign();
503 den = zero;
504 return *this;
505 }
506
507 // Protect against self-modification
508 IntType r_num = r.num;
509 IntType r_den = r.den;
510
511 // Avoid overflow and preserve normalization
514 num = (num/gcd1) * (r_num/gcd2);
515 den = (den/gcd2) * (r_den/gcd1);
516 return *this;
517}
518
519template <typename IntType>
521{
522 IntType zero(0);
523
524 // handle the Inf and NaN cases
525 if(den == zero)
526 {
527 if(r.den == zero)
528 throw bad_rational();
529 if(r.num != zero)
530 num *= r.sign();
531 return *this;
532 }
533 if(r.num == zero)
534 {
535 if(num == zero)
536 throw bad_rational();
537 num = IntType(sign()); // normalized inf!
538 den = zero;
539 return *this;
540 }
541
542 if (num == zero)
543 return *this;
544
545 // Protect against self-modification
546 IntType r_num = r.num;
547 IntType r_den = r.den;
548
549 // Avoid overflow and preserve normalization
552 num = (num/gcd1) * (r_den/gcd2);
553 den = (den/gcd2) * (r_num/gcd1);
554
555 if (den < zero) {
556 num = -num;
557 den = -den;
558 }
559 return *this;
560}
561
562// Mixed-mode operators -- implement explicitly to save gcd() calculations
563template <typename IntType>
564inline Rational<IntType>&
566{
567 num += i * den;
568 return *this;
569}
570
571template <typename IntType>
572inline Rational<IntType>&
574{
575 num -= i * den;
576 return *this;
577}
578
579template <typename IntType>
582{
583 if(i == IntType(1))
584 return *this;
585 IntType zero(0);
586 if(i == zero)
587 {
588 if(den == zero)
589 {
590 throw bad_rational();
591 }
592 else
593 {
594 num = zero;
595 den = IntType(1);
596 return *this;
597 }
598 }
599
600 IntType g = gcd(i, den);
601 den /= g;
602 num *= i / g;
603 return *this;
604}
605
606template <typename IntType>
609{
610 if(i == IntType(1))
611 return *this;
612
613 IntType zero(0);
614 if(i == zero)
615 {
616 if(num == zero)
617 throw bad_rational();
618 num = IntType(sign()); // normalized inf!
619 den = zero;
620 return *this;
621 }
622
623 IntType g = gcd(i, num);
624 if(i < zero)
625 {
626 num /= -g;
627 den *= -i / g;
628 }
629 else
630 {
631 num /= g;
632 den *= i / g;
633 }
634 return *this;
635}
636
637// Increment and decrement
638template <typename IntType>
640{
641 // This can never denormalise the fraction
642 num += den;
643 return *this;
644}
645
646template <typename IntType>
648{
649 // This can never denormalise the fraction
650 num -= den;
651 return *this;
652}
653
654// Normalisation
655template <typename IntType>
657{
658 // Avoid repeated construction
659 IntType zero(0);
660
661 if (den == zero)
662 {
663 if(num == zero)
664 throw bad_rational();
665 if(num < zero)
666 num = IntType(-1);
667 else
668 num = IntType(1);
669 return;
670 }
671
672 // Handle the case of zero separately, to avoid division by zero
673 if (num == zero) {
674 den = IntType(1);
675 return;
676 }
677
678 IntType g = gcd<IntType>(num, den);
679
680 num /= g;
681 den /= g;
682
683 // Ensure that the denominator is positive
684 if (den < zero) {
685 num = -num;
686 den = -den;
687 }
688}
689
690/********************************************************/
691/* */
692/* Rational-Traits */
693/* */
694/********************************************************/
695
696/** \page RationalTraits Numeric and Promote Traits of Rational
697
698 The numeric and promote traits for Rational follow
699 the general specifications for \ref NumericPromotionTraits and
700 \ref AlgebraicField. They are implemented in terms of the traits of the basic types by
701 partial template specialization:
702
703 \code
704
705 template <class T>
706 struct NumericTraits<Rational<T> >
707 {
708 typedef Rational<typename NumericTraits<T>::Promote> Promote;
709 typedef Rational<typename NumericTraits<T>::RealPromote> RealPromote;
710
711 typedef typename NumericTraits<T>::isIntegral isIntegral;
712 typedef VigraTrueType isScalar;
713 typedef typename NumericTraits<T>::isSigned isSigned;
714 typedef VigraTrueType isOrdered;
715
716 // etc.
717 };
718
719 template<class T>
720 struct NormTraits<Rational<T> >
721 {
722 typedef Rational<T> Type;
723 typedef typename NumericTraits<Type>::Promote SquaredNormType;
724 typedef Type NormType;
725 };
726
727 template <class T1, class T2>
728 struct PromoteTraits<Rational<T1>, Rational<T2> >
729 {
730 typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
731 };
732 \endcode
733
734 <b>\#include</b> <vigra/rational.hxx><br>
735 Namespace: vigra
736
737*/
738#ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
739
740template<class T>
741struct NumericTraits<Rational<T> >
742{
743 typedef Rational<T> Type;
744 typedef Rational<typename NumericTraits<T>::Promote> Promote;
745 typedef Rational<typename NumericTraits<T>::RealPromote> RealPromote;
746 typedef std::complex<Rational<RealPromote> > ComplexPromote;
747 typedef T ValueType;
748
749 typedef typename NumericTraits<T>::isIntegral isIntegral;
750 typedef VigraTrueType isScalar;
751 typedef typename NumericTraits<T>::isSigned isSigned;
752 typedef VigraTrueType isOrdered;
753 typedef VigraFalseType isComplex;
754
755 static Type zero() { return Type(0); }
756 static Type one() { return Type(1); }
757 static Type nonZero() { return one(); }
758
759 static Promote toPromote(Type const & v)
760 { return Promote(v.numerator(), v.denominator(), false); }
761 static RealPromote toRealPromote(Type const & v)
762 { return RealPromote(v.numerator(), v.denominator(), false); }
763 static Type fromPromote(Promote const & v)
764 { return Type(NumericTraits<T>::fromPromote(v.numerator()),
765 NumericTraits<T>::fromPromote(v.denominator()), false); }
766 static Type fromRealPromote(RealPromote const & v)
767 { return Type(NumericTraits<T>::fromRealPromote(v.numerator()),
768 NumericTraits<T>::fromRealPromote(v.denominator()), false); }
769};
770
771template<class T>
772struct NormTraits<Rational<T> >
773{
774 typedef Rational<T> Type;
775 typedef typename NumericTraits<Type>::Promote SquaredNormType;
776 typedef Type NormType;
777};
778
779template <class T>
780struct PromoteTraits<Rational<T>, Rational<T> >
781{
782 typedef Rational<typename PromoteTraits<T, T>::Promote> Promote;
783 static Promote toPromote(Rational<T> const & v) { return v; }
784};
785
786template <class T1, class T2>
787struct PromoteTraits<Rational<T1>, Rational<T2> >
788{
789 typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
790 static Promote toPromote(Rational<T1> const & v) { return v; }
791 static Promote toPromote(Rational<T2> const & v) { return v; }
792};
793
794template <class T1, class T2>
795struct PromoteTraits<Rational<T1>, T2 >
796{
797 typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
798 static Promote toPromote(Rational<T1> const & v) { return v; }
799 static Promote toPromote(T2 const & v) { return Promote(v); }
800};
801
802template <class T1, class T2>
803struct PromoteTraits<T1, Rational<T2> >
804{
805 typedef Rational<typename PromoteTraits<T1, T2>::Promote> Promote;
806 static Promote toPromote(T1 const & v) { return Promote(v); }
807 static Promote toPromote(Rational<T2> const & v) { return v; }
808};
809
810#endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
811
812/********************************************************/
813/* */
814/* RationalOperations */
815/* */
816/********************************************************/
817
818/** \addtogroup RationalOperations Functions for Rational
819
820 \brief <b>\#include</b> <vigra/rational.hxx><br>
821
822 These functions fulfill the requirements of an \ref AlgebraicField.
823
824 Namespace: vigra
825 <p>
826
827 */
828//@{
829
830/********************************************************/
831/* */
832/* arithmetic */
833/* */
834/********************************************************/
835
836 /// unary plus
837template <typename IntType>
839{
840 return r;
841}
842
843 /// unary minus (negation)
844template <typename IntType>
846{
847 return Rational<IntType>(-r.numerator(), r.denominator(), false);
848}
849
850 /// addition
851template <typename IntType>
853{
854 return l += r;
855}
856
857 /// subtraction
858template <typename IntType>
860{
861 return l -= r;
862}
863
864 /// multiplication
865template <typename IntType>
867{
868 return l *= r;
869}
870
871 /// division
872template <typename IntType>
874{
875 return l /= r;
876}
877
878 /// addition of right-hand <tt>IntType</tt> argument
879template <typename IntType>
880inline Rational<IntType>
882{
883 return l += r;
884}
885
886 /// subtraction of right-hand <tt>IntType</tt> argument
887template <typename IntType>
888inline Rational<IntType>
890{
891 return l -= r;
892}
893
894 /// multiplication with right-hand <tt>IntType</tt> argument
895template <typename IntType>
896inline Rational<IntType>
898{
899 return l *= r;
900}
901
902 /// division by right-hand <tt>IntType</tt> argument
903template <typename IntType>
904inline Rational<IntType>
906{
907 return l /= r;
908}
909
910 /// addition of left-hand <tt>IntType</tt> argument
911template <typename IntType>
912inline Rational<IntType>
914{
915 return r += l;
916}
917
918 /// subtraction from left-hand <tt>IntType</tt> argument
919template <typename IntType>
920inline Rational<IntType>
922{
923 return (-r) += l;
924}
925
926 /// multiplication with left-hand <tt>IntType</tt> argument
927template <typename IntType>
928inline Rational<IntType>
930{
931 return r *= l;
932}
933
934 /// division of left-hand <tt>IntType</tt> argument
935template <typename IntType>
936inline Rational<IntType>
938{
939 if(r.numerator() < IntType(0))
940 return Rational<IntType>(-r.denominator(), -r.numerator(), false) *= l;
941 else
942 return Rational<IntType>(r.denominator(), r.numerator(), false) *= l;
943}
944
945/********************************************************/
946/* */
947/* comparison */
948/* */
949/********************************************************/
950
951
952 /// equality
953template <typename IntType1, typename IntType2>
954inline bool
956{
957 return l.denominator() == r.denominator() &&
958 l.numerator() == r.numerator(); // works since numbers are normalized, even
959 // if they represent +-infinity
960}
961
962 /// equality with right-hand <tt>IntType2</tt> argument
963template <typename IntType1, typename IntType2>
964inline bool
966{
967 return ((l.denominator() == IntType1(1)) && (l.numerator() == i));
968}
969
970 /// equality with left-hand <tt>IntType1</tt> argument
971template <typename IntType1, typename IntType2>
972inline bool
974{
975 return r == l;
976}
977
978 /// inequality
979template <typename IntType1, typename IntType2>
980inline bool
982{
983 return l.denominator() != r.denominator() ||
984 l.numerator() != r.numerator(); // works since numbers are normalized, even
985 // if they represent +-infinity
986}
987
988 /// inequality with right-hand <tt>IntType2</tt> argument
989template <typename IntType1, typename IntType2>
990inline bool
992{
993 return ((l.denominator() != IntType1(1)) || (l.numerator() != i));
994}
995
996 /// inequality with left-hand <tt>IntType1</tt> argument
997template <typename IntType1, typename IntType2>
998inline bool
1000{
1001 return r != l;
1002}
1003
1004 /// less-than
1005template <typename IntType1, typename IntType2>
1006bool
1008{
1009 // Avoid repeated construction
1010 typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
1011 IntType zero(0);
1012
1013 // Handle the easy cases. Take advantage of the fact
1014 // that the denominator is never negative.
1015 if(l.denominator() == zero)
1016 {
1017 if(r.denominator() == zero)
1018 // -inf < inf, !(-inf < -inf), !(inf < -inf), !(inf < inf)
1019 return l.numerator() < r.numerator();
1020 else
1021 // -inf < -1, -inf < 0, -inf < 1
1022 // !(inf < -1), !(inf < 0), !(inf < 1)
1023 return l.numerator() < zero;
1024 }
1025 if(r.denominator() == zero)
1026 // -1 < inf, 0 < inf, 1 < inf
1027 // !(-1 < -inf), !(0 < -inf), !(1 < -inf)
1028 return r.numerator() > zero;
1029 // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
1030 if(l.numerator() >= zero && r.numerator() <= zero)
1031 return false;
1032 // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
1033 if(l.numerator() <= zero && r.numerator() >= zero)
1034 return true;
1035
1036 // both numbers have the same sign (and are neither zero or +-infinity)
1037 // => calculate result, avoid overflow
1038 IntType gcd1 = gcd<IntType>(l.numerator(), r.numerator());
1039 IntType gcd2 = gcd<IntType>(r.denominator(), l.denominator());
1040 return (l.numerator()/gcd1) * (r.denominator()/gcd2) <
1041 (l.denominator()/gcd2) * (r.numerator()/gcd1);
1042}
1043
1044 /// less-than with right-hand <tt>IntType2</tt> argument
1045template <typename IntType1, typename IntType2>
1046bool
1047operator< (const Rational<IntType1> & l, IntType2 const & i)
1048{
1049 // Avoid repeated construction
1050 typedef typename PromoteTraits<IntType1, IntType2>::Promote IntType;
1051 IntType zero(0);
1052
1053 // Handle the easy cases. Take advantage of the fact
1054 // that the denominator is never negative.
1055 if(l.denominator() == zero)
1056 // -inf < -1, -inf < 0, -inf < 1
1057 // !(inf < -1), !(inf < 0), !(inf < 1)
1058 return l.numerator() < zero;
1059 // !(1 < -1), !(1 < 0), !(0 < -1), !(0 < 0)
1060 if(l.numerator() >= zero && i <= zero)
1061 return false;
1062 // -1 < 0, -1 < 1, 0 < 1 (note: !(0 < 0) was already handled!)
1063 if(l.numerator() <= zero && i >= zero)
1064 return true;
1065
1066 // Now, use the fact that n/d truncates towards zero as long as n and d
1067 // are both positive.
1068 // Divide instead of multiplying to avoid overflow issues. Of course,
1069 // division may be slower, but accuracy is more important than speed...
1070 if (l.numerator() > zero)
1071 return (l.numerator()/l.denominator()) < i;
1072 else
1073 return -i < (-l.numerator()/l.denominator());
1074}
1075
1076 /// less-than with left-hand <tt>IntType1</tt> argument
1077template <typename IntType1, typename IntType2>
1078inline bool
1079operator<(IntType1 const & l, Rational<IntType2> const & r)
1080{
1081 return r > l;
1082}
1083
1084 /// greater-than
1085template <typename IntType1, typename IntType2>
1086inline bool
1088{
1089 return r < l;
1090}
1091
1092 /// greater-than with right-hand <tt>IntType2</tt> argument
1093template <typename IntType1, typename IntType2>
1094bool
1096{
1097 // Trap equality first
1098 if (l.numerator() == i && l.denominator() == IntType1(1))
1099 return false;
1100
1101 // Otherwise, we can use operator<
1102 return !(l < i);
1103}
1104
1105 /// greater-than with left-hand <tt>IntType1</tt> argument
1106template <typename IntType1, typename IntType2>
1107inline bool
1109{
1110 return r < l;
1111}
1112
1113 /// less-equal
1114template <typename IntType1, typename IntType2>
1115inline bool
1117{
1118 return !(r < l);
1119}
1120
1121 /// less-equal with right-hand <tt>IntType2</tt> argument
1122template <typename IntType1, typename IntType2>
1123inline bool
1124operator<=(Rational<IntType1> const & l, IntType2 const & r)
1125{
1126 return !(l > r);
1127}
1128
1129 /// less-equal with left-hand <tt>IntType1</tt> argument
1130template <typename IntType1, typename IntType2>
1131inline bool
1132operator<=(IntType1 const & l, Rational<IntType2> const & r)
1133{
1134 return r >= l;
1135}
1136
1137 /// greater-equal
1138template <typename IntType1, typename IntType2>
1139inline bool
1141{
1142 return !(l < r);
1143}
1144
1145 /// greater-equal with right-hand <tt>IntType2</tt> argument
1146template <typename IntType1, typename IntType2>
1147inline bool
1149{
1150 return !(l < r);
1151}
1152
1153 /// greater-equal with left-hand <tt>IntType1</tt> argument
1154template <typename IntType1, typename IntType2>
1155inline bool
1157{
1158 return r <= l;
1159}
1160
1161/********************************************************/
1162/* */
1163/* algebraic functions */
1164/* */
1165/********************************************************/
1166
1167 /// absolute value
1168template <typename IntType>
1169inline Rational<IntType>
1171{
1172 if (r.numerator() >= IntType(0))
1173 return r;
1174
1175 return Rational<IntType>(-r.numerator(), r.denominator(), false);
1176}
1177
1178 /// norm (same as <tt>abs(r)</tt>)
1179template <typename IntType>
1180inline Rational<IntType>
1182{
1183 return abs(r);
1184}
1185
1186 /// squared norm
1187template <typename IntType>
1188inline typename NormTraits<Rational<IntType> >::SquaredNormType
1190{
1191 return typename NormTraits<Rational<IntType> >::SquaredNormType(sq(r.numerator()), sq(r.denominator()), false);
1192}
1193
1194 /** integer powers
1195
1196 <tt>throws bad_rational</tt> if indeterminate expression.
1197 */
1198template <typename IntType>
1199Rational<IntType>
1200pow(const Rational<IntType>& r, int e)
1201{
1202 IntType zero(0);
1203 int ae;
1204 if(e == 0)
1205 {
1206 if(r.denominator() == zero)
1207 throw bad_rational();
1208 return Rational<IntType>(IntType(1));
1209 }
1210 else if(e < 0)
1211 {
1212 if(r.numerator() == zero)
1213 return Rational<IntType>(IntType(1), zero, false);
1214 if(r.denominator() == zero)
1215 return Rational<IntType>(zero);
1216 ae = -e;
1217 }
1218 else
1219 {
1220 if(r.denominator() == zero || r.numerator() == zero)
1221 return r;
1222 ae = e;
1223 }
1224
1225 IntType nold = r.numerator(), dold = r.denominator(),
1226 nnew = IntType(1), dnew = IntType(1);
1227 for(; ae != 0; ae >>= 1, nold *= nold, dold *= dold)
1228 {
1229 if(ae % 2 != 0)
1230 {
1231 nnew *= nold;
1232 dnew *= dold;
1233 }
1234 }
1235 if(e < 0)
1236 {
1237 if(nnew < zero)
1238 return Rational<IntType>(-dnew, -nnew, false);
1239 else
1240 return Rational<IntType>(dnew, nnew, false);
1241 }
1242 else
1243 return Rational<IntType>(nnew, dnew, false);
1244}
1245
1246 /// largest integer not larger than <tt>r</tt>
1247template <typename IntType>
1248Rational<IntType>
1250{
1251 IntType zero(0), one(1);
1252 if(r.denominator() == zero || r.denominator() == one)
1253 return r;
1254 return r.numerator() < zero ?
1255 Rational<IntType>(r.numerator() / r.denominator() - one)
1256 : Rational<IntType>(r.numerator() / r.denominator());
1257}
1258
1259 /// smallest integer not smaller than <tt>r</tt>
1260template <typename IntType>
1261Rational<IntType>
1263{
1264 IntType zero(0), one(1);
1265 if(r.denominator() == zero || r.denominator() == one)
1266 return r;
1267 return r.numerator() < IntType(0) ?
1268 Rational<IntType>(r.numerator() / r.denominator())
1269 : Rational<IntType>(r.numerator() / r.denominator() + one);
1270}
1271
1272
1273 /** Type conversion
1274
1275 Executes <tt>static_cast<T>(numerator()) / denominator()</tt>.
1276
1277 <b>Usage:</b>
1278
1279 \code
1280 Rational<int> r;
1281 int i;
1282 double d;
1283 i = rational_cast<int>(r); // round r downwards
1284 d = rational_cast<double>(r); // represent rational as a double
1285 r = rational_cast<Rational<int> >(r); // no change
1286 \endcode
1287 */
1288template <typename T, typename IntType>
1290{
1291 return static_cast<T>(src.numerator())/src.denominator();
1292}
1293
1294template <class T>
1295inline T const & rational_cast(T const & v)
1296{
1297 return v;
1298}
1299
1300//@}
1301
1302template <typename IntType>
1303std::ostream& operator<< (std::ostream& os, const vigra::Rational<IntType>& r)
1304{
1305 os << r.numerator() << '/' << r.denominator();
1306 return os;
1307}
1308
1309} // namespace vigra
1310
1311#endif // VIGRA_RATIONAL_HPP
1312
Class for a single RGB value.
Definition rgbvalue.hxx:128
RGBValue()
Definition rgbvalue.hxx:209
Definition rational.hxx:194
Rational(param_type n, param_type d, bool doNormalize=true)
Definition rational.hxx:236
Rational(param_type n)
Definition rational.hxx:223
param_type numerator() const
Definition rational.hxx:269
Rational & operator=(param_type n)
Definition rational.hxx:261
Rational(double v, double epsilon=1e-4)
Definition rational.hxx:248
bool is_ninf() const
Definition rational.hxx:351
bool operator!() const
Definition rational.hxx:339
Rational & operator/=(const Rational &r)
Definition rational.hxx:520
param_type denominator() const
Definition rational.hxx:273
bool is_pinf() const
Definition rational.hxx:343
Rational operator++(int)
Definition rational.hxx:332
Rational & operator--()
Definition rational.hxx:647
bool is_inf() const
Definition rational.hxx:359
IntType value_type
Definition rational.hxx:198
Rational & operator+=(const Rational &r)
Definition rational.hxx:402
Rational & operator*=(const Rational &r)
Definition rational.hxx:486
Rational operator--(int)
Definition rational.hxx:335
Rational(Rational< U > const &r)
Definition rational.hxx:216
Rational()
Definition rational.hxx:208
If< typenameTypeTraits< IntType >::isBuiltinType, IntType, IntTypeconst & >::type param_type
Definition rational.hxx:204
Rational & operator++()
Definition rational.hxx:639
Rational & assign(param_type n, param_type d, bool doNormalize=true)
Definition rational.hxx:391
Rational & operator-=(const Rational &r)
Definition rational.hxx:452
int sign() const
Definition rational.hxx:369
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition fixedpoint.hxx:667
Diff2D operator+(Diff2D const &a, Diff2D const &b)
Definition diff2d.hxx:725
bool operator<=(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
less or equal
Definition fixedpoint.hxx:521
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
bool operator>=(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater or equal
Definition fixedpoint.hxx:539
IntType lcm(IntType n, IntType m)
Definition rational.hxx:122
IntType gcd(IntType n, IntType m)
Definition rational.hxx:81
T sign(T t)
The sign function.
Definition mathutil.hxx:591
bool operator>(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater
Definition fixedpoint.hxx:530
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
bool operator<(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
less than
Definition fixedpoint.hxx:512
T rational_cast(const Rational< IntType > &src)
Definition rational.hxx:1289
bool operator==(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
equal
Definition fftw3.hxx:825
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition fixedpoint.hxx:675
bool operator!=(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
not equal
Definition fftw3.hxx:841
Diff2D operator-(Diff2D const &a, Diff2D const &b)
Definition diff2d.hxx:697

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.12.1