9dde43a4dad9006713ccd819c265adb200d886b3
[ipdf/code.git] / rational.h
1 #ifndef _RATIONAL_H
2 #define _RATIONAL_H
3
4 /**
5  * A really shoddy implementation of Rational numbers
6  */
7
8 #include "common.h"
9 #include <cmath>
10 #include <cassert>
11
12 namespace IPDF
13 {
14
15 /* Recursive version  of GCD
16 template <class T>
17 T gcd(const T & a, const T & b)
18 {
19         Debug("Called on %li/%li", int64_t(a), int64_t(b));
20         if (a == T(1) || a == T(0)) return T(1);
21         if (b == T(0)) return a;
22         if (b == a) 
23         {
24                 Debug("Equal!");
25                 return a;
26         }
27         Debug("Not equal!");
28         
29         if (a > b) return gcd(a-b,b);
30         return gcd(a, b-a);
31 }
32 */
33
34 /** Greatest Common Divisor of p and q **/
35 template <class T>
36 T gcd(const T & p, const T & q)
37 {
38
39
40         T g(1);
41         T big(p);
42         T small(q);
43         if (p < q)
44         {
45                 big = q;
46                 small = p;
47         }
48         if (small == T(0))
49                 return g;
50         while ((g = big % small) > T(0))
51         {
52                 //Debug("big = %li, small = %li", int64_t(big), int64_t(small));
53                 big = small;
54                 small = g;
55                 //Debug("Loop %u", ++count);
56         }
57         return small;
58 }       
59
60
61 template <class T = int64_t>
62 struct Rational
63 {
64         /** Construct from a double.**/
65         Rational(double d=0) : P(d*1e6), Q(1e6) // Possibly the worst thing ever...
66         {
67                 Simplify();
68                 if (!CheckAccuracy(d, "Construct from double"))
69                 {
70                         //Fatal("Bwah bwah :(");
71                 }
72         }
73
74         Rational(const T & _P, const T & _Q) : P(_P), Q(_Q)
75         {
76                 Simplify();
77         }
78
79         Rational(const Rational & cpy) : P(cpy.P), Q(cpy.Q)
80         {
81                 Simplify();
82         }
83
84         void Simplify()
85         {
86                 if (Q < T(0)) 
87                 {
88                         P = -P;
89                         Q = -Q;
90                 }
91                 if (P == T(0))
92                 {
93                         Q = T(1);
94                         return;
95                 }
96                 T g = gcd(T(llabs(P)),T(llabs(Q)));
97                 //Debug("Got gcd!");
98                 P /= g;
99                 Q /= g;
100         }
101
102         bool operator==(const Rational & r)  const
103         {
104                 if (P == r.P && Q == r.Q) return true;
105                 return ToDouble() == r.ToDouble();
106         }
107
108
109         bool operator<(const Rational & r) const {return (P*r.Q < r.P * Q);}
110         bool operator>(const Rational & r) const {return !(*this < r);}
111         bool operator<=(const Rational & r) const {return *this == r || *this < r;}
112         bool operator>=(const Rational & r) const {return *this == r || *this > r;}
113         bool operator!=(const Rational & r) const {return !(*this == r);}
114
115         Rational operator+(const Rational & r) const 
116         {
117                 Rational result = (r.P == T(0)) ? Rational(P,Q) : Rational(P*r.Q + r.P*Q, Q*r.Q);
118                 if (!result.CheckAccuracy(ToDouble() * r.ToDouble(),"+"))
119                 {
120                         Debug("This is %s (%f) and r is %s (%f)", Str().c_str(), ToDouble(), r.Str().c_str(), r.ToDouble());
121                 }
122                 return result;
123         }
124         Rational operator-(const Rational & r) const 
125         {
126                 Rational result = (r.P == T(0)) ? Rational(P,Q) : Rational(P*r.Q - r.P*Q, Q*r.Q);
127                 result.CheckAccuracy(ToDouble() - r.ToDouble(),"-");
128                 return result;
129         }
130         Rational operator*(const Rational & r) const 
131         {
132                 Rational result(P * r.P, Q * r.Q);
133                 if (!result.CheckAccuracy(ToDouble() * r.ToDouble(),"*"))
134                 {
135                         Debug("This is %s (%f) and r is %s (%f)", Str().c_str(), ToDouble(), r.Str().c_str(), r.ToDouble());
136                 }
137                 return result;
138         }
139         Rational operator/(const Rational & r) const 
140         {
141                 Rational result(P * r.Q, Q*r.P);
142                 if (!result.CheckAccuracy(ToDouble() / r.ToDouble(),"/"))
143                 {
144                         Debug("This is %s (%f) and r is %s (%f)", Str().c_str(), ToDouble(), r.Str().c_str(), r.ToDouble());
145                 }
146                 return result;
147         }       
148
149         /** To cheat, use these **/
150         //Rational operator+(const Rational & r) const {return Rational(ToDouble()+r.ToDouble());}
151         //Rational operator-(const Rational & r) const {return Rational(ToDouble()-r.ToDouble());}
152         //Rational operator*(const Rational & r) const {return Rational(ToDouble()*r.ToDouble());}
153         //Rational operator/(const Rational & r) const {return Rational(ToDouble()/r.ToDouble());}
154
155         Rational & operator=(const Rational & r) {P = r.P; Q = r.Q; return *this;}
156         Rational & operator+=(const Rational & r) {this->operator=(*this+r); return *this;}
157         Rational & operator-=(const Rational & r) {this->operator=(*this-r); return *this;}
158         Rational & operator*=(const Rational & r) {this->operator=(*this*r); return *this;}
159         Rational & operator/=(const Rational & r) {this->operator=(*this/r); return *this;}
160
161         double ToDouble() const {return (double)(P) / (double)(Q);}
162         bool CheckAccuracy(double d, const char * msg, double threshold = 1e-3) const
163         {
164                 double result = fabs(ToDouble() - d);
165                 if (d != 0e0) result /= d;
166                 if (result > threshold)
167                 {
168                         Warn("(%s) : Rational %s (%f) is not close enough at representing %f (%f vs %f)", msg, Str().c_str(), ToDouble(), d, result, threshold);
169                         Backtrace();
170                         return false;
171                 }
172                 return true;
173         }
174         std::string Str() const
175         {
176                 std::stringstream s;
177                 s << int64_t(P) << "/" << int64_t(Q);
178                 return s.str();
179         }
180         
181         T P;
182         T Q;
183 };
184
185 inline Rational<int64_t> pow(const Rational<int64_t> & a, const Rational<int64_t> & b)
186 {
187         //TODO:Implement properly
188         int64_t P = std::pow((double)a.P, b.ToDouble());
189         int64_t Q = std::pow((double)a.Q, b.ToDouble());
190         return Rational<int64_t>(P, Q);
191 }
192
193
194
195 }
196
197 #endif //_RATIONAL_H
198

UCC git Repository :: git.ucc.asn.au