00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_RATIONAL_HPP
00014 #define __MMX_RATIONAL_HPP
00015 #include <numerix/integer.hpp>
00016 namespace mmx {
00017
00018 class rational;
00019 class rational_rep;
00020 rational copy (const rational& i);
00021
00022
00023
00024
00025
00026 class rational_rep REP_STRUCT {
00027 private:
00028 mpq_t x;
00029 public:
00030 inline rational_rep () { mpq_init (x); }
00031 inline ~rational_rep () { mpq_clear (x); }
00032 friend class rational;
00033 template<typename V> friend class floating;
00034 };
00035
00036 class rational {
00037 INDIRECT_PROTO (rational, rational_rep)
00038 public:
00039 inline const mpq_t& operator * () const { return rep->x; }
00040 inline mpq_t& operator * () { secure(); return rep->x; }
00041
00042 public:
00043
00044 inline rational (): rep (new rational_rep ()) {}
00045 inline rational (const format<rational>& fm): rep (new rational_rep ()) {
00046 (void) fm; }
00047 inline rational (signed int i): rep (new rational_rep ()) {
00048 mpq_set_si (rep->x, i, 1); }
00049 inline rational (unsigned int i):
00050 rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
00051 inline rational (signed short int i):
00052 rep (new rational_rep ()) { mpq_set_si (rep->x, i, 1); }
00053 inline rational (unsigned short int i):
00054 rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
00055 inline rational (signed long int i):
00056 rep (new rational_rep ()) { mpq_set_si (rep->x, i, 1); }
00057 inline rational (unsigned long int i):
00058 rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
00059 inline rational (const integer& i):
00060 rep (new rational_rep ()) { mpq_set_z (rep->x, *i); }
00061 inline rational (const double& d):
00062 rep (new rational_rep ()) { mpq_set_d (rep->x, d); }
00063
00064 inline rational (char * const & c):
00065 rep (new rational_rep ()) { mpq_set_str (rep->x, c, 10); }
00066
00067
00068 inline friend integer numerator (const rational& x1) {
00069 integer r; mpq_get_num (*r, *x1); return r; }
00070 inline friend integer denominator (const rational& x1) {
00071 integer r; mpq_get_den (*r, *x1); return r; }
00072 inline friend bool is_integer (const rational& x) {
00073 return mpz_cmp_ui (mpq_denref (*x), 1) == 0; }
00074
00075
00076 inline friend rational copy (const rational& x1) {
00077 rational r; mpq_set (*r, *x1); return r; }
00078 inline friend rational operator - (const rational& x1) {
00079 rational r; mpq_neg (*r, *x1); return r; }
00080 inline friend rational operator + (const rational& x1, const rational& x2) {
00081 rational r; mpq_add (*r, *x1, *x2); return r; }
00082 inline friend rational operator - (const rational& x1, const rational& x2) {
00083 rational r; mpq_sub (*r, *x1, *x2); return r; }
00084 inline friend rational operator * (const rational& x1, const rational& x2) {
00085 rational r; mpq_mul (*r, *x1, *x2); return r; }
00086 inline friend rational square (const rational& x1) {
00087 rational r; mpq_mul (*r, *x1, *x1); return r; }
00088 inline friend rational operator / (const rational& x1, const rational& x2) {
00089 ASSERT (mpq_sgn (*x2) != 0, "division by zero");
00090 rational r; mpq_div (*r, *x1, *x2); return r; }
00091 inline friend bool divides (const rational& x1, const rational& x2) {
00092 (void) x2;
00093 return mpq_sgn (*x1) != 0; }
00094 inline friend bool is_invertible (const rational& x) {
00095 return x != 0; }
00096 inline friend rational invert (const rational& x) {
00097 return rational (1) / x; }
00098 inline friend rational operator << (const rational& x1, const xint& shift) {
00099 rational r; mpq_mul_2si (*r, *x1, shift); return r; }
00100 inline friend rational operator >> (const rational& x1, const xint& shift) {
00101 rational r; mpq_mul_2si (*r, *x1, -shift); return r; }
00102 inline friend rational abs (const rational& x1) {
00103 rational r; mpq_abs (*r, *x1); return r; }
00104
00105
00106 inline friend rational& operator += (rational& x1, const rational& x2) {
00107 x1.secure (); mpq_add (*x1, *x1, *x2); return x1; }
00108 inline friend rational& operator -= (rational& x1, const rational& x2) {
00109 x1.secure (); mpq_sub (*x1, *x1, *x2); return x1; }
00110 inline friend rational& operator *= (rational& x1, const rational& x2) {
00111 x1.secure (); mpq_mul (*x1, *x1, *x2); return x1; }
00112 inline friend rational& operator /= (rational& x1, const rational& x2) {
00113 ASSERT (mpq_sgn (*x2) != 0, "division by zero");
00114 x1.secure (); mpq_div (*x1, *x1, *x2); return x1; }
00115 inline friend rational& operator <<= (rational& x1, const xint& shift) {
00116 x1.secure (); mpq_mul_2si (*x1, *x1, shift); return x1; }
00117 inline friend rational& operator >>= (rational& x1, const xint& shift) {
00118 x1.secure (); mpq_mul_2si (*x1, *x1, -shift); return x1; }
00119
00120 inline friend void
00121 neg (rational& r) {
00122 r.secure (); mpq_neg (*r, *r); }
00123 inline friend void
00124 neg (rational& r, const rational& x1) {
00125 r.secure (); mpq_neg (*r, *x1); }
00126 inline friend void
00127 add (rational& r, const rational& x1, const rational& x2) {
00128 r.secure (); mpq_add (*r, *x1, *x2); }
00129 inline friend void
00130 sub (rational& r, const rational& x1, const rational& x2) {
00131 r.secure (); mpq_sub (*r, *x1, *x2); }
00132 inline friend void
00133 mul (rational& r, const rational& x1, const rational& x2) {
00134 r.secure (); mpq_mul (*r, *x1, *x2); }
00135 inline friend void
00136 div (rational& r, const rational& x1, const rational& x2) {
00137 ASSERT (mpq_sgn (*x2) != 0, "division by zero");
00138 r.secure (); mpq_div (*r, *x1, *x2); }
00139
00140
00141 inline friend bool operator == (const rational& x1, const rational& x2) {
00142 return mpq_cmp (*x1, *x2) == 0; }
00143 inline friend bool operator != (const rational& x1, const rational& x2) {
00144 return mpq_cmp (*x1, *x2) != 0; }
00145 inline friend bool operator < (const rational& x1, const rational& x2) {
00146 return mpq_cmp (*x1, *x2) < 0; }
00147 inline friend bool operator > (const rational& x1, const rational& x2) {
00148 return mpq_cmp (*x1, *x2) > 0; }
00149 inline friend bool operator <= (const rational& x1, const rational& x2) {
00150 return mpq_cmp (*x1, *x2) <= 0; }
00151 inline friend bool operator >= (const rational& x1, const rational& x2) {
00152 return mpq_cmp (*x1, *x2) >= 0; }
00153 inline friend int sign (const rational& x1) {
00154 int r= mpq_sgn (*x1); return (r<0? -1: (r==0? 0: 1)); }
00155 inline friend rational min (const rational& x1, const rational& x2) {
00156 return x1 <= x2? x1: x2; }
00157 inline friend rational max (const rational& x1, const rational& x2) {
00158 return x1 >= x2? x1: x2; }
00159
00160
00161 inline friend rational floor (const rational& x1) {
00162 return quo (numerator (x1), denominator (x1)); }
00163 inline friend rational trunc (const rational& x1) {
00164 return sign (x1) * quo (abs (numerator (x1)), denominator (x1)); }
00165 inline friend rational ceil (const rational& x1) {
00166 return quo (numerator (x1) + denominator (x1) - 1, denominator (x1)); }
00167 inline friend rational round (const rational& x1) {
00168 return quo (numerator (x1) + (denominator (x1) >> 1), denominator (x1)); }
00169
00170
00171 inline friend bool is_square (const rational& a) {
00172 return mpz_perfect_square_p (mpq_numref (*a)) &&
00173 mpz_perfect_square_p (mpq_denref (*a)); }
00174 inline friend rational sqrt (const rational& a) {
00175 VERIFY (is_square (a), "not a perfect square");
00176 rational r;
00177 mpz_sqrt (mpq_numref (*r), mpq_numref (*a));
00178 mpz_sqrt (mpq_denref (*r), mpq_denref (*a));
00179 return r; }
00180 inline friend rational pow (const rational& a, const int& i) {
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 if (i >= 0)
00195 return rational (pow (numerator (a), i)) / pow (denominator (a), i);
00196 else
00197 return rational (pow (denominator (a), -i)) / pow (numerator (a), -i);
00198 }
00199 inline friend rational pow (const rational& a, const integer& i) {
00200 ASSERT (is_int (i) || a == rational (0) ||
00201 a == rational (1) || a == rational (-1), "power too large");
00202 return pow (a, as_int (i)); }
00203
00204
00205 inline friend nat hash (const rational& x) {
00206 const __mpq_struct &rep= (*x)[0];
00207 if (rep._mp_num._mp_size == 0) return 0;
00208 return
00209 ((nat) (rep._mp_num._mp_d[0])) ^
00210 ((nat) (rep._mp_num._mp_size << 3)) ^
00211 ((nat) (rep._mp_den._mp_d[0] << 7)) ^
00212 ((nat) (rep._mp_den._mp_size << 11));
00213 }
00214
00215 template<typename V> friend class floating;
00216 template<typename V> friend floating<V> as_floating (const rational& x);
00217 friend double as_double (const rational& q) {return mpq_get_d(*q);}
00218 };
00219 INDIRECT_IMPL (rational, rational_rep)
00220
00221 template<> struct symbolic_type_information<rational> {
00222 static const nat id= SYMBOLIC_RATIONAL; };
00223
00224 UNARY_RETURN_TYPE(STMPL,numerator_op,rational,integer);
00225 UNARY_RETURN_TYPE(STMPL,denominator_op,rational,integer);
00226
00227
00228
00229
00230
00231 syntactic flatten (const rational& x);
00232
00233 template<>
00234 struct binary_helper<rational>: public void_binary_helper<rational> {
00235 static inline string short_type_name () { return "Q"; }
00236 static inline generic full_type_name () { return "Rational"; }
00237 static inline generic disassemble (const rational& x) {
00238 return gen_vec (as<generic> (numerator (x)),
00239 as<generic> (denominator (x))); }
00240 static inline rational assemble (const generic& v) {
00241 return rational (as<integer> (vector_access (v, 0))) /
00242 rational (as<integer> (vector_access (v, 1))); }
00243 static inline void write (const port& out, const rational& x) {
00244 binary_write<integer> (out, numerator (x));
00245 binary_write<integer> (out, denominator (x)); }
00246 static inline rational read (const port& in) {
00247 integer n= binary_read<integer> (in);
00248 integer d= binary_read<integer> (in);
00249 return rational (n) / rational (d); }
00250 };
00251
00252 template<>
00253 struct as_helper<double,rational> {
00254 static inline double cv (const rational& x) { return as_double (x); }
00255 };
00256
00257 template<typename F> inline rational
00258 promote (const F& x, const rational&) { return as<rational> (x); }
00259
00260
00261
00262
00263
00264 generic construct (const rational& x);
00265 template<> inline rational duplicate (const rational& x) { return copy (x); }
00266 inline void set_default (rational& r) { r= 0; }
00267 inline void set_nan (rational& r) { r= 0; }
00268
00269 TRUE_TO_EXACT_IDENTITY_SUGAR(,rational)
00270 ARITH_INT_SUGAR(,rational)
00271 GCD_SUGAR(,rational)
00272 POOR_MAN_ELEMENTARY_SUGAR(,rational)
00273
00274
00275
00276
00277
00278
00279 #define rational_new(n,d) rational(n)/d
00280 #define old_integer_pow(i,j) (j>=0?as<generic>(pow(i,j)):as<generic>(1/rational(pow(i,-j))))
00281 inline rational integer_pow (const integer& a, const int& i) {
00282 if (i >= 0) return pow (a, i); else return rational (1) / pow (a, -i); }
00283 inline rational integer_pow (const integer& a, const integer& i) {
00284 if (i >= 0) return pow (a, i); else return rational (1) / pow (a, -i); }
00285
00286 }
00287 #endif // __MMX_RATIONAL_HPP