/Users/mourrain/Devel/mmx/numerix/include/numerix/rational.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : rational.hpp
* DESCRIPTION: Rational numbers
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_RATIONAL_HPP
#define __MMX_RATIONAL_HPP
#include <numerix/integer.hpp>
namespace mmx {

class rational;
class rational_rep;
rational copy (const rational& i);

/******************************************************************************
* Rational numbers
******************************************************************************/

class rational_rep REP_STRUCT {
private:
  mpq_t x;
public:
  inline rational_rep () { mpq_init (x); }
  inline ~rational_rep () { mpq_clear (x); }
  friend class rational;
  template<typename V> friend class floating;
};

class rational {
INDIRECT_PROTO (rational, rational_rep)
private:
  inline const mpq_t& operator * () const { return rep->x; }
  inline mpq_t& operator * () { secure(); return rep->x; }

public:
  // constructors
  inline rational (): rep (new rational_rep ()) {}
  inline rational (signed int i): rep (new rational_rep ()) {
    mpq_set_si (rep->x, i, 1); }
  inline rational (unsigned int i):
    rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
  inline rational (signed short int i):
    rep (new rational_rep ()) { mpq_set_si (rep->x, i, 1); }
  inline rational (unsigned short int i):
    rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
  inline rational (signed long int i):
    rep (new rational_rep ()) { mpq_set_si (rep->x, i, 1); }
  inline rational (unsigned long int i):
    rep (new rational_rep ()) { mpq_set_ui (rep->x, i, 1); }
  inline rational (const integer& i):
    rep (new rational_rep ()) { mpq_set_z (rep->x, *i); }
  inline rational (const double& d):
    rep (new rational_rep ()) { mpq_set_d (rep->x, d); }

  inline rational (char * const & c):
    rep (new rational_rep ()) { mpq_set_str (rep->x, c, 10); }

  // getting the numerator and denominator
  inline friend integer numerator (const rational& x1) {
    integer r; mpq_get_num (*r, *x1); return r; }
  inline friend integer denominator (const rational& x1) {
    integer r; mpq_get_den (*r, *x1); return r; }

  // basic arithmetic
  inline friend rational copy (const rational& x1) {
    rational r; mpq_set (*r, *x1); return r; }
  inline friend rational operator - (const rational& x1) {
    rational r; mpq_neg (*r, *x1); return r; }
  inline friend rational operator + (const rational& x1, const rational& x2) {
    rational r; mpq_add (*r, *x1, *x2); return r; }
  inline friend rational operator - (const rational& x1, const rational& x2) {
    rational r; mpq_sub (*r, *x1, *x2); return r; }
  inline friend rational operator * (const rational& x1, const rational& x2) {
    rational r; mpq_mul (*r, *x1, *x2); return r; }
  inline friend rational square (const rational& x1) {
    rational r; mpq_mul (*r, *x1, *x1); return r; }
  inline friend rational operator / (const rational& x1, const rational& x2) {
    ASSERT (mpq_sgn (*x2) != 0, "division by zero");
    rational r; mpq_div (*r, *x1, *x2); return r; }
  inline friend bool divides (const rational& x1, const rational& x2) {
    return mpq_sgn (*x1) != 0; }
  inline friend rational invert (const rational& x) {
    return rational (1) / x; }
  inline friend rational operator << (const rational& x1, const xint& shift) {
    rational r; mpq_mul_2si (*r, *x1, shift); return r; }
  inline friend rational operator >> (const rational& x1, const xint& shift) {
    rational r; mpq_mul_2si (*r, *x1, -shift); return r; }
  inline friend rational abs (const rational& x1) {
    rational r; mpq_abs (*r, *x1); return r; }

  // in-place arithmetic
  inline friend rational& operator += (rational& x1, const rational& x2) {
    x1.secure (); mpq_add (*x1, *x1, *x2); return x1; }
  inline friend rational& operator -= (rational& x1, const rational& x2) {
    x1.secure (); mpq_sub (*x1, *x1, *x2); return x1; }
  inline friend rational& operator *= (rational& x1, const rational& x2) {
    x1.secure (); mpq_mul (*x1, *x1, *x2); return x1; }
  inline friend rational& operator /= (rational& x1, const rational& x2) {
    ASSERT (mpq_sgn (*x2) != 0, "division by zero");
    x1.secure (); mpq_div (*x1, *x1, *x2); return x1; }
  inline friend rational& operator <<= (rational& x1, const xint& shift) {
    x1.secure (); mpq_mul_2si (*x1, *x1, shift); return x1; }
  inline friend rational& operator >>= (rational& x1, const xint& shift) {
    x1.secure (); mpq_mul_2si (*x1, *x1, -shift); return x1; }

  inline friend void
  neg (rational& r) {
    r.secure (); mpq_neg (*r, *r); }
  inline friend void
  neg (rational& r, const rational& x1) {
    r.secure (); mpq_neg (*r, *x1); }
  inline friend void
  add (rational& r, const rational& x1, const rational& x2) {
    r.secure (); mpq_add (*r, *x1, *x2); }
  inline friend void
  sub (rational& r, const rational& x1, const rational& x2) {
    r.secure (); mpq_sub (*r, *x1, *x2); }
  inline friend void
  mul (rational& r, const rational& x1, const rational& x2) {
    r.secure (); mpq_mul (*r, *x1, *x2); }
  inline friend void
  div (rational& r, const rational& x1, const rational& x2) {
    ASSERT (mpq_sgn (*x2) != 0, "division by zero");
    r.secure (); mpq_div (*r, *x1, *x2); }

  // comparing rational numbers
  inline friend bool operator == (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) == 0; }
  inline friend bool operator != (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) != 0; }
  inline friend bool operator < (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) < 0; }
  inline friend bool operator > (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) > 0; }
  inline friend bool operator <= (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) <= 0; }
  inline friend bool operator >= (const rational& x1, const rational& x2) {
    return mpq_cmp (*x1, *x2) >= 0; }
  inline friend int sign (const rational& x1) {
    int r= mpq_sgn (*x1); return (r<0? -1: (r==0? 0: 1)); }
  inline friend rational min (const rational& x1, const rational& x2) {
    return x1 <= x2? x1: x2; }
  inline friend rational max (const rational& x1, const rational& x2) {
    return x1 >= x2? x1: x2; }

  // rounding towards integers
  inline friend rational floor (const rational& x1) {
    return quo (numerator (x1), denominator (x1)); }
  inline friend rational trunc (const rational& x1) {
    return sign (x1) * quo (abs (numerator (x1)), denominator (x1)); }
  inline friend rational ceil (const rational& x1) {
    return quo (numerator (x1) + denominator (x1) - 1, denominator (x1)); }
  inline friend rational round (const rational& x1) {
    return quo (numerator (x1) + (denominator (x1) >> 1), denominator (x1)); }

  // other functions on rationals
  inline friend bool is_square (const rational& a) {
    return mpz_perfect_square_p (mpq_numref (*a)) &&
           mpz_perfect_square_p (mpq_denref (*a)); }
  inline friend rational sqrt (const rational& a) {
    VERIFY (is_square (a), "not a perfect square");
    rational r;
    mpz_sqrt (mpq_numref (*r), mpq_numref (*a));
    mpz_sqrt (mpq_denref (*r), mpq_denref (*a));
    return r; }
  inline friend rational pow (const rational& a, const int& i) {
    rational r;
    if (i >= 0) {
      mpz_pow_ui (mpq_numref (*r), mpq_numref (*a), i);
      mpz_pow_ui (mpq_denref (*r), mpq_denref (*a), i);
    }
    else {
      ASSERT (mpq_sgn (*a) != 0, "division by zero");
      mpz_pow_ui (mpq_denref (*r), mpq_numref (*a), -i);
      mpz_pow_ui (mpq_numref (*r), mpq_denref (*a), -i);
    }
    return r; }
  inline friend rational pow (const rational& a, const integer& i) {
    // FIXME: test whether integer i fits in int
    return pow (a, as_int (i)); }

  // miscellaneous functions
  inline friend nat hash (const rational& x) {
    const __mpq_struct &rep= (*x)[0];
    if (rep._mp_num._mp_size == 0) return 0;
    return
      ((nat) (rep._mp_num._mp_d[0])) ^
      ((nat) (rep._mp_num._mp_size << 3)) ^
      ((nat) (rep._mp_den._mp_d[0] << 7)) ^
      ((nat) (rep._mp_den._mp_size << 11));
  }

  template<typename V> friend class floating;
  template<typename V> friend floating<V> as_floating (const rational& x);
  friend double as_double (const rational& q)  {return mpq_get_d(*q);}
};
INDIRECT_IMPL (rational, rational_rep)

template<> struct symbolic_type_information<rational> {
  static const nat id= SYMBOLIC_RATIONAL; };

UNARY_RETURN_TYPE(STMPL,numerator_op,rational,integer);
UNARY_RETURN_TYPE(STMPL,denominator_op,rational,integer);

/******************************************************************************
* Conversion related
******************************************************************************/

syntactic flatten (const rational& x);

template<>
struct binary_helper<rational>: public void_binary_helper<rational> {
  static inline string short_type_name () { return "Q"; }
  static inline generic full_type_name () { return "Rational"; }
  static inline generic disassemble (const rational& x) {
    return gen_vec (as<generic> (numerator (x)),
                    as<generic> (denominator (x))); }
  static inline rational assemble (const generic& v) {
    return rational (as<integer> (vector_access (v, 0))) /
           rational (as<integer> (vector_access (v, 1))); }
  static inline void write (const port& out, const rational& x) {
    binary_write<integer> (out, numerator (x));
    binary_write<integer> (out, denominator (x)); }
  static inline rational read (const port& in) {
    integer n= binary_read<integer> (in);
    integer d= binary_read<integer> (in);
    return rational (n) / rational (d); }
};

template<>
struct as_helper<double,rational> {
  static inline double cv (const rational& x) { return as_double (x); }
};

template<typename F> inline rational
promote (const F& x, const rational&) { return as<rational> (x); }

/******************************************************************************
* Defaults
******************************************************************************/

generic construct (const rational& x);
template<> inline rational duplicate (const rational& x) { return copy (x); }
inline void set_nan (rational& r) { r= 0; }

TRUE_TO_EXACT_IDENTITY_SUGAR(,rational)
ARITH_INT_SUGAR(,rational)
GCD_SUGAR(,rational)
POOR_MAN_ELEMENTARY_SUGAR(,rational)

/******************************************************************************
* Extras for glue
******************************************************************************/

#define rational_new(n,d) rational(n)/d
#define integer_pow(i,j) (j>=0?as<generic>(pow(i,j)):as<generic>(1/rational(pow(i,-j))))
  // FIXME: to be removed

}
#endif // __MMX_RATIONAL_HPP