atlas  0.6
arithmetic.h
Go to the documentation of this file.
1 // This is arithmetic.h
2 /*
3  Copyright (C) 2004,2005 Fokko du Cloux
4  Copyright (C) 2006-2016 Marc van Leeuwen
5  part of the Atlas of Lie Groups and Representations
6
8 */
9
10 #ifndef ARITHMETIC_H /* guard against multiple inclusions */
11 #define ARITHMETIC_H
12
13 #include "arithmetic_fwd.h"
14
15 #include <iostream>
16
17 /******** function declarations **********************************************/
18
19 namespace atlas {
20
21 namespace arithmetic {
22
23 // operations without conversion for all integral types (all are inlined below)
24  template<typename I> I divide(I, Denom_t);
25  template<typename I> Denom_t remainder(I, Denom_t);
26  // template<typename I> I factorial(I); // moved to the size module
27
29
30  Denom_t gcd (Numer_t, Denom_t); // signed first argument only!
31
32  // the following functions are entirely unsigned
34
35  Denom_t div_gcd (Denom_t a, Denom_t b); // $a/\gcd(a,b)$
36
37  Denom_t lcm (Denom_t a, Denom_t b,
38  Denom_t& gcd=dummy_gcd, Denom_t& a_mult=dummy_mult);
39  // after |lcm(a,b,d,p)| one has |d=gcd(a,b)| and |p%a==0|, |p%b==d|, |p<lcm|
40
41  Denom_t power(Denom_t base, unsigned int exponent);
42
43  // inlined; first two arguments are supposed already reduced modulo third
45
47
48 } // |namespace arithmetic|
49
50 /* Class definitions */
51
52 namespace arithmetic {
53
54 class Rational
55 {
58 public:
59  explicit Rational (Numer_t n=0,Denom_t d=1) : num(n), denom(d)
60  { normalize(); }
61
62  Numer_t numerator() const { return num; }
63
64  /* the C++ rule that one unsigned operand silently converts the other
65  operand to unsigned as well makes exporting denominator as unsigned too
66  error prone; e.g., floor=numerator()/denominator() would wreak havoc */
67  Numer_t denominator() const { return Numer_t(denom); }
68
69  // these operators all return normalised results
70  Rational operator+(Rational q) const;
71  Rational operator-(Rational q) const;
72  Rational operator*(Rational q) const;
73  Rational operator/(Rational q) const; // assumes $q\neq0$, will not throw
74  Rational operator%(Rational q) const; // assumes $q\neq0$, will not throw
75
77  { num=q.num; denom=q.denom; return normalize(); }
78
79  Rational& operator+=(Rational q) { return operator=(operator+(q)); }
80  Rational& operator-=(Rational q) { return operator=(operator-(q)); }
81  Rational& operator*=(Rational q) { return operator=(operator*(q)); }
82  Rational& operator/=(Rational q) { return operator=(operator/(q)); }
83  Rational& operator%=(Rational q) { return operator=(operator%(q)); }
84
85  // assignment operators with integers have efficient implemementations
88  Rational& operator*=(Numer_t n);
89  Rational& operator/=(Numer_t n); // assumes $n\neq0$, will not throw
90  Rational& operator%=(Numer_t n); // assumes $n\neq0$, will not throw
91
92  // these definitions must use |denominator()| to ensure signed comparison
93  bool operator==(Rational q) const
94  { return num*q.denominator()==denominator()*q.num; }
95  bool operator!=(Rational q) const
96  { return num*q.denominator()!=denominator()*q.num; }
97  bool operator<(Rational q) const
98  { return num*q.denominator()<denominator()*q.num; }
99  bool operator<=(Rational q) const
100  { return num*q.denominator()<=denominator()*q.num; }
101  bool operator>(Rational q) const
102  { return num*q.denominator()>denominator()*q.num; }
103  bool operator>=(Rational q) const
104  { return num*q.denominator()>=denominator()*q.num; }
105
106  inline Rational& normalize();
107  Rational& power(int n); // raise to power |n| and return |*this|
108
109 }; // |class Rational|
110
111
113 {
114  int real_part, s_part;
115  public:
116  explicit Split_integer(int a=0, int b=0) : real_part(a), s_part(b) {}
117
118  int e() const { return real_part; }
119  int s() const { return s_part; }
120
121  bool operator== (Split_integer y) const { return e()==y.e() and s()==y.s(); }
122  bool operator!= (Split_integer y) const { return not operator==(y); }
123
125  { real_part+=y.e(); s_part+=y.s(); return *this; }
127  { real_part-=y.e(); s_part-=y.s(); return *this; }
128  Split_integer operator +(Split_integer y) { return y+= *this; }
129  Split_integer operator -(Split_integer y) { return y.negate()+= *this; }
130  Split_integer operator -() const { return Split_integer(*this).negate(); }
132  { return Split_integer(e()*y.e()+s()*y.s(),e()*y.s()+s()*y.e()); }
133
134  Split_integer& operator*= (int n) { real_part*=n; s_part*=n; return *this; }
135  Split_integer& operator*= (Split_integer y) { return *this=operator*(y); }
136  Split_integer& negate() { real_part=-e(); s_part=-s(); return *this; }
137  Split_integer& times_s() { std::swap(real_part,s_part); return *this; }
138  Split_integer& times_1_s() { real_part= -(s_part-=e()); return *this; }
139
140 }; // |class Split_integer|
141
142 std::ostream& operator<< (std::ostream& out, const Rational& frac);
143
144
145 /******** inline function definitions ***************************************/
146
163 template<typename I>
164  inline I divide(I a, Denom_t b)
165  { return a >= 0 ? a/b : ~(~a/b); } // left operand is safely made unsigned
166
167 // override for |I=Denom_t| (is instantiated for |Matrix<Denom_t>|)
169  { return a/b; } // unsigned division is OK in this case
170
171 /*
172  Return the remainder of the division of |a| (signed) by |b| (unsigned).
173
174  The point is to allow |I| to be a signed type, and avoid the catastrophic
175  implicit conversion to unsigned when using the '%' operation. Also corrects
176  the deficiency of 'signed modulo' by always returning the unique number |r|
177  in [0,m[ such that $a = q.b + r$, in other words with |q=divide(a,b)| above.
178
179  NOTE: For $a<0$ one should \emph{not} return |b - (-a % b)|; this fails when
180  $b$ divides $a$. However replacing |-| by |~|, which maps $a\mapsto-1-a$
181  and satifies |~(q*b+r)==~q*b+(b+~r)|, the result is always correct.
182 */
183 template<typename I>
184  inline Denom_t remainder(I a, Denom_t b)
185  { return a >= 0 ? a%b // safe implicit conversion to unsigned here
186  : b+~(~static_cast<Denom_t>(a)%b); // safe explicit conversion here
187  }
188
189  inline Denom_t div_gcd (Denom_t d, Denom_t a) { return d/unsigned_gcd(a,d); }
190
191  inline Denom_t gcd (Numer_t a, Denom_t b)
192  {
193  if (a < 0)
194  return unsigned_gcd(static_cast<Denom_t>(-a),b);
195  else
196  return unsigned_gcd(static_cast<Denom_t>(a),b);
197  }
198
199  // we assume |a| and |b| to be less than |n| here
201  {
202  if (a < n-b)
203  return a + b;
204  else
205  return a - (n-b);
206  }
207
209  {
210  Denom_t d = gcd(num,denom);
211  if (d>1) num/=Numer_t(d),denom/=d;
212  return *this;
213  }
214
215 } // |namespace arithmetic|
216
217 } // |namespace atlas|
218
219 #endif
Denom_t dummy_mult
Definition: arithmetic.cpp:50
Denom_t remainder(I, Denom_t)
Definition: arithmetic.h:184
Rational & operator-=(Rational q)
Definition: arithmetic.h:80
bool operator>(Rational q) const
Definition: arithmetic.h:101
Vector< C > operator/(Vector< C > v, C c)
Definition: matrix.h:136
bool operator!=(const type_expr &x, const type_expr &y)
Definition: axis-types.h:374
Denom_t gcd(Numer_t, Denom_t)
Definition: arithmetic.h:191
bool operator<(Rational q) const
Definition: arithmetic.h:97
Rational & operator%=(Rational q)
Definition: arithmetic.h:83
Rational & operator=(Rational q)
Definition: arithmetic.h:76
PID_Matrix< C > operator-(PID_Matrix< C > A, C c)
Definition: matrix.h:51
Definition: arithmetic.h:54
bool operator==(Rational q) const
Definition: arithmetic.h:93
Denom_t modProd(Denom_t a, Denom_t b, Denom_t n)
Definition: arithmetic.cpp:85
Denom_t power(Denom_t x, unsigned int n)
Definition: arithmetic.cpp:93
Numer_t num
Definition: arithmetic.h:56
Denom_t div_gcd(Denom_t a, Denom_t b)
Definition: arithmetic.h:189
I divide(I, Denom_t)
Definition: arithmetic.h:164
Split_integer & times_s()
Definition: arithmetic.h:137
Vector< C > & operator%=(Vector< C > &v, C c)
Definition: matrix.cpp:127
void swap(simple_list< T, Alloc > &x, simple_list< T, Alloc > &y)
Definition: sl_list.h:674
Numer_t numerator() const
Definition: arithmetic.h:62
Rational & operator/=(Rational q)
Definition: arithmetic.h:82
PID_Matrix< C > operator+(PID_Matrix< C > A, C c)
Definition: matrix.h:44
Split_integer & negate()
Definition: arithmetic.h:136
Rational(Numer_t n=0, Denom_t d=1)
Definition: arithmetic.h:59
bool operator<=(Rational q) const
Definition: arithmetic.h:99
RationalVector< C2 > operator*(const matrix::Matrix< C1 > &M, const RationalVector< C2 > &v)
Definition: ratvec.cpp:158
Split_integer & times_1_s()
Definition: arithmetic.h:138
PID_Matrix< C > & operator+=(PID_Matrix< C > &A, C c)
Definition: matrix.cpp:733
#define out(c)
Definition: cweave.c:205
int s() const
Definition: arithmetic.h:119
Denom_t denom
Definition: arithmetic.h:57
std::ostream & operator<<(std::ostream &strm, const Split_integer &s)
Definition: basic_io.cpp:157
PID_Matrix< C > & operator-=(PID_Matrix< C > &A, C c)
Definition: matrix.h:48
int s_part
Definition: arithmetic.h:114
unsigned long long int Denom_t
Definition: Atlas.h:69
Vector< C > & operator/=(Vector< C > &v, C c)
Definition: matrix.cpp:99
Split_integer(int a=0, int b=0)
Definition: arithmetic.h:116
long long int Numer_t
Definition: Atlas.h:68
Rational & operator+=(Rational q)
Definition: arithmetic.h:79
Denom_t dummy_gcd
Definition: arithmetic.cpp:50
Rational & normalize()
Definition: arithmetic.h:208
Denom_t unsigned_gcd(Denom_t a, Denom_t b)
Definition: arithmetic.cpp:40
Permutation normalize(const DynkinDiagram &d)
Definition: dynkin.cpp:386
unsigned long n
Definition: axis.cpp:77
bool operator==(const type_expr &x, const type_expr &y)
Definition: axis-types.cpp:257
bool operator>=(Rational q) const
Definition: arithmetic.h:103
Definition: Atlas.h:38
Definition: arithmetic.h:200
int e() const
Definition: arithmetic.h:118
const expr & e
Definition: axis.cpp:95
Definition: arithmetic.h:112
Denom_t lcm(Denom_t a, Denom_t b, Denom_t &gcd, Denom_t &mult_a)
Definition: arithmetic.cpp:58
Rational & operator*=(Rational q)
Definition: arithmetic.h:81
bool operator!=(Rational q) const
Definition: arithmetic.h:95
Numer_t denominator() const
Definition: arithmetic.h:67