include/algebramix/quotient.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : quotient.hpp
* DESCRIPTION: Fractions
* COPYRIGHT  : (C) 2007  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* 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_QUOTIENT__HPP
#define __MMX_QUOTIENT__HPP
#include <basix/port.hpp>
#include <basix/function.hpp>
namespace mmx {
struct std_quotient {};
#define TMPL_DEF template<typename NT, typename DT>
#define TMPL template<typename NT, typename DT>
#define TMPLK template<typename NT, typename DT, typename K>
#define Quotient quotient<NT,DT>
TMPL class quotient;
TMPL inline NT numerator (const Quotient& x);
TMPL inline DT denominator (const Quotient& x);

/******************************************************************************
* Normalization helper
******************************************************************************/

TMPL_DEF
struct quotient_normalization_helper {
  static inline void normalize (NT& a, DT& b) {}
};

template<typename C,typename V> struct polynomial;

template<typename C,typename V>
struct quotient_normalization_helper<polynomial<C,V>,
                                     polynomial<C,V> > {
  static inline void normalize (polynomial<C,V>& a,
                                polynomial<C,V>& b) {
    ASSERT (b != 0, "unexpected null denominator");
    C c= b [N(b)-1];
    if (is_invertible (c)) { a= a / c; b= b / c; } }
};

/******************************************************************************
* The quotient type
******************************************************************************/

TMPL_DEF
class quotient {
MMX_ALLOCATORS;
private:
  NT n;
  DT d;
public:
  inline quotient ():
    n (), d (1) {}
  template<typename T> inline quotient (const format<T>& fm):
    n (promote (0, fm)), d (promote (1, fm)) {}
  template<typename T> inline quotient (const T& x):
    n (outplace_set_as<NT> (x)), d (outplace_set_as<DT> ((int) 1)) {}
  template<typename NT2, typename DT2>
  inline quotient (const quotient<NT2,DT2>& x):
    n (numerator (x)), d (denominator (x)) {}
  template<typename NT2, typename DT2> inline
  quotient (const NT2& x, const DT2& y):
    n (x), d (y) {
      if (n == promote (0, n)) d= promote (1, d);
      else quotient_normalization_helper<NT,DT>::normalize (n, d); }
  template<typename NT2, typename DT2> inline
  quotient (const NT2& x, const DT2& y, bool simplify):
    n (x), d (y) {
    assert (simplify);
    if (n == promote (0, n))
      d= promote (1, d);
    else {
      DT g= gcd (n, d);
      n /= g;
      d /= g;
    }
    quotient_normalization_helper<NT,DT>::normalize (n, d); }
  friend NT numerator LESSGTR (const Quotient& x);
  friend DT denominator LESSGTR (const Quotient& x);
};
DEFINE_UNARY_FORMAT_2 (quotient)

STYPE_TO_TYPE(TMPL,scalar_type,Quotient,NT);
UNARY_RETURN_TYPE(TMPL,numerator_op,Quotient,NT);
UNARY_RETURN_TYPE(TMPL,denominator_op,Quotient,DT);

/******************************************************************************
* Basic routines
******************************************************************************/

TMPL inline NT numerator   (const Quotient& x) { return x.n; }
TMPL inline DT denominator (const Quotient& x) { return x.d; }

TMPL inline format<NT> CF (const Quotient& x) {
  return get_format (numerator (x)); }

TMPL inline nat precision (const Quotient& x) {
  return max (precision (numerator (x)), precision (denominator (x))); }

/* FIXME: why does this not work? A compiler bug in gcc 2.95.3?
template<typename Op, typename NT, typename DT> inline nat
unary_hash (const Quotient& x) {
  nat h= Op::op (numerator (x));
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ Op::op (denominator (x));
}

template<typename Op, typename NT, typename DT> bool
binary_test (const Quotient& x1, const Quotient& x2) {
  return Op::op (numerator (x1), numerator (x2)) &&
         Op::op (denominator (x1), denominator (x2));
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Quotient)
EXACT_IDENTITY_OP_SUGAR(TMPL,Quotient)
HARD_IDENTITY_OP_SUGAR(TMPL,Quotient)
*/

template<typename S1, typename D> quotient<D,D>
map (const function_1<D,Argument(S1) >& fun, const quotient<S1,S1>& x) {
  return quotient<D,D> (fun (numerator (x)), fun (denominator (x)));
}

template<typename S1, typename D> quotient<D,D>
map (D (*fun) (const S1&), const quotient<S1,S1>& x) {
  return quotient<D,D> (fun (numerator (x)), fun (denominator (x)));
}

TMPL inline nat hash (const Quotient& x) {
  nat h= hash (numerator (x));
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ hash (denominator (x)); }
TMPL inline nat exact_hash (const Quotient& x) {
  nat h= exact_hash (numerator (x));
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ exact_hash (denominator (x)); }
TMPL inline nat hard_hash (const Quotient& x) {
  nat h= hard_hash (numerator (x));
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ hard_hash (denominator (x)); }
TMPL inline bool operator == (const Quotient& x1, const Quotient& x2) {
  return numerator (x1) * denominator (x2) ==
         numerator (x2) * denominator (x1); }
TMPL inline bool operator != (const Quotient& x1, const Quotient& x2) {
  return !(x1 == x2); }
TMPL inline bool exact_eq (const Quotient& x1, const Quotient& x2) {
  return exact_eq (numerator (x1), numerator (x2)) &&
         exact_eq (denominator (x1), denominator (x2)); }
TMPL inline bool exact_neq (const Quotient& x1, const Quotient& x2) {
  return !exact_eq (x1, x2); }
TMPL inline bool hard_eq (const Quotient& x1, const Quotient& x2) {
  return hard_eq (numerator (x1), numerator (x2)) &&
         hard_eq (denominator (x1), denominator (x2)); }
TMPL inline bool hard_neq (const Quotient& x1, const Quotient& x2) {
  return !hard_eq (x1, x2); }

EQUAL_INT_SUGAR(TMPL,Quotient);
EQUAL_SCALAR_SUGAR(TMPL,Quotient,NT)
EQUAL_SCALAR_SUGAR_BIS(TMPL,Quotient,NT)

TMPL syntactic
flatten (const Quotient& x) {
  return flatten (numerator (x)) / flatten (denominator (x));
}

TMPL
struct binary_helper<Quotient >: public void_binary_helper<Quotient > {
  static inline string short_type_name () {
    return "Qu" * Short_type_name (NT) * Short_type_name (DT); }
  static inline generic full_type_name () {
    return gen ("Quotient", Full_type_name (NT), Full_type_name (DT)); }
  static inline generic disassemble (const Quotient& q) {
    return gen_vec (as<generic> (numerator (q)),
                    as<generic> (denominator (q))); }
  static inline Quotient assemble (const generic& v) {
    return Quotient (as<NT> (vector_access (v, 0)),
                     as<DT> (vector_access (v, 1))); }
  static inline void write (const port& out, const Quotient& q) {
    binary_write<NT> (out, numerator (q));
    binary_write<DT> (out, denominator (q)); }
  static inline Quotient read (const port& in) {
    NT n= binary_read<NT> (in);
    DT d= binary_read<DT> (in);
    return Quotient (n, d); }
};

GCD_SUGAR(TMPL,Quotient)

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Quotient
operator + (const Quotient& x1, const Quotient& x2) {
  DT g  = gcd (denominator (x1), denominator(x2));
  DT y1 = denominator (x1) / g;
  DT y2 = denominator (x2) / g;
  NT n  = numerator(x1) * y2 + numerator(x2) * y1;
  DT h  = gcd (g, n);
  return Quotient (n / h, (g / h) * y1 * y2);
}

TMPL inline Quotient
operator + (const Quotient& x1, const NT& x2) {
  return Quotient (numerator (x1) + x2 * denominator (x1),
                   denominator (x1));
}

TMPL inline Quotient
operator + (const NT& x1, const Quotient& x2) {
  return Quotient (x1 * denominator (x2) + numerator (x2),
                   denominator (x2));
}

TMPL inline Quotient
operator - (const Quotient& x) {
  return Quotient (-numerator (x), denominator (x));
}

TMPL inline Quotient
operator - (const Quotient& x1, const Quotient& x2) {
  DT g  = gcd (denominator (x1), denominator(x2));
  DT y1 = denominator (x1) / g;
  DT y2 = denominator (x2) / g;
  NT n  = numerator(x1) * y2 - numerator(x2) * y1;
  DT h  = gcd (g, n);
  return Quotient (n / h, (g / h) * y1 * y2);
}

TMPL inline Quotient
operator - (const Quotient& x1, const NT& x2) {
  return Quotient (numerator (x1) - x2 * denominator (x1),
                   denominator (x1));
}

TMPL inline Quotient
operator - (const NT& x1, const Quotient& x2) {
  return Quotient (x1 * denominator (x2) - numerator (x2),
                   denominator (x2));
}

TMPL inline Quotient
operator * (const NT& c, const Quotient& x) {
  DT g = gcd (c, denominator (x));
  return Quotient ((c / g) * numerator (x), denominator (x) / g);
}

TMPL inline Quotient
operator * (const Quotient& x, const NT& c) {
  DT g = gcd (c, denominator (x));
  return Quotient (numerator (x) * (c / g), denominator (x) / g);
}

TMPL inline Quotient
operator * (const Quotient& x1, const Quotient& x2) {
  DT g1 = gcd (numerator (x1), denominator (x2));
  DT g2 = gcd (numerator (x2), denominator (x1));
  return Quotient ((numerator (x1) / g1) * (numerator (x2) / g2),
                   (denominator (x1) / g2) * (denominator (x2) / g1));
}

TMPL inline Quotient
operator / (const NT& c, const Quotient& x) {
  // assumes NT = DT
  ASSERT (numerator (x) != 0, "division by zero");
  DT g= gcd (c, numerator (x));
  return Quotient ((c / g) * denominator (x), numerator (x) / g);
}

TMPL inline Quotient
operator / (const Quotient& x, const DT& c) {
  ASSERT (c != 0, "division by zero");
  DT g= gcd (numerator (x), c);
  return Quotient (numerator (x) / g, denominator (x) * (c / g));
}

TMPL inline Quotient
operator / (const Quotient& x1, const Quotient& x2) {
  // assumes NT = DT
  ASSERT (numerator (x2) != 0, "division by zero");
  NT g1= gcd (numerator (x1), numerator (x2));
  DT g2= gcd (denominator (x1), denominator (x2));
  return Quotient ((numerator (x1) / g1) * (denominator (x2) / g2),
                   (denominator (x1) / g2) * (numerator (x2) / g1));
}

TMPL inline bool
is_invertible (const Quotient& x) {
  return !is_exact_zero (x);  // assumes NT = DT
}

ARITH_SCALAR_INT_SUGAR(TMPL,Quotient);

/******************************************************************************
* Ordering
******************************************************************************/

TMPL inline int
sign (const Quotient& x) {
  return sign (numerator (x)) * sign (denominator (x));
}

COMPARE_SUGAR(TMPL,Quotient);
COMPARE_INT_SUGAR(TMPL,Quotient);
COMPARE_SCALAR_SUGAR(TMPL,Quotient,NT)
COMPARE_SCALAR_SUGAR_BIS(TMPL,Quotient,NT)

/******************************************************************************
* Other operations
******************************************************************************/

TMPL Quotient
derive (const Quotient& x) {
  return Quotient (derive (numerator (x)) * denominator (x) -
                   numerator (x) * derive (denominator (x)),
                   square (denominator (x)),
                   true);
}

TMPL Quotient
xderive (const Quotient& x) {
  return Quotient (xderive (numerator (x)) * denominator (x) -
                   numerator (x) * xderive (denominator (x)),
                   square (denominator (x)),
                   true);
}

template<typename NT, typename DT, typename VT> Quotient
derive (const Quotient& x, const VT& v) {
  return Quotient (derive (numerator (x), v) * denominator (x) -
                   numerator (x) * derive (denominator (x), v),
                   square (denominator (x)),
                   true);
}

/******************************************************************************
* Other operations
******************************************************************************/

BINARY_RETURN_TYPE(TMPLK,evaluate_op,Quotient,K,
                   Binary_return_type(evaluate_op,NT,K));

TMPLK Binary_return_type(evaluate_op,Quotient,K)
evaluate (const Quotient& x, const K& a) {
  return evaluate (numerator (x), a) / evaluate (denominator (x), a); } 

TMPLK bool
is_evaluable (const Quotient& x, const K& a, Evaluate_type(Quotient,K)& b) {
  Evaluate_type(NT,K) num= b;
  Evaluate_type(DT,K) den= b;
  if (!is_evaluable (numerator (x), a, num)) return false;
  if (!is_evaluable (denominator (x), a, den)) return false;
  if (is_exact_zero (den)) return false; // FIXME: test for invertibility
  b= num / den; return true; }

#undef TMPL_DEF
#undef TMPL
#undef TMPLK
#undef Quotient
} // namespace mmx
#endif // __MMX_QUOTIENT__HPP

---