Exact Arithmetic using GMP

Prerequisits

#include <Integer.h>
#include <Rational.h>

Most calculations in polymake are made exactly, using rational numbers. The best implementation known to us which provides arbitrary precision and high-tuned performance on various computer platforms, is GMP. On the other hand, we did like very much the convenient interface of the Rational and Integer classes from the (in the meanwhile obsolete and no more supported) libg++, which allowed to use the rational values in expressions in the most natural way, as if they were built-in numeric types. So we have decided to union the advantages of both and have written thin wrapper classes mimicking the old Rational and Integer. They perform almost no calculations on their own, but delegate the whole hard job to GMP functions.

In the meanwhile, GMP has got its own C++ wrapper classes. They are implemented differently from our classes: all arithmetic operators are "lazy", returning expression templates instead of ready results. While this technique has proven to improve the performance in longer chained expressions, and is massively used in polymake vector and matrix classes too, we consider it a bit dangerous for the basic numerical type.

In an innocent expression like (a+b)*M, where a and b are rational scalars and M is a rational matrix, the sum a+b would be repeatedly calculated for each element of the resulting matrix. Since such mixed expressions are prevailing in polymake code, we have decided to stay with our own wrappers, relying on the so called return value optimization of the compiler. Future changes in GMP may let us revise this decision.

Constructors, conversions

GMP numbers can be constructed from the most built-in types. They can be converted to each other and to the most built-in types.

One of Integer, int, or long.

Any combination of Integer, int, and long.

Any combination of Rational, Integer, int, and long.

The real return type is a temporary proxy object forwarding all access operations to the internal representation of the numerator rsp. denominator and taking care of the data consistency.

Integer();

Rational();

Initialize with 0.

Integer(INT i);

Rational(INT i);

Integer::operator=(INT i);

Rational::operator=(INT i);

Initialize with an integral number.

Integer(double d);

Integer::operator=(double d);

Integer(const Rational& r);

Integer::operator=(double d);

Initialize with the integral part of a floating-point number without rounding.

Rational(double d);

Rational(INT num, INT den) throw(gmp_error);

Rational::set(INT num, INT den) throw(gmp_error);

Initialize with a fractional number. An attempt to assign zero to the denominator raises an exception.

Integer::operator int() const throw(gmp_error);

Integer::operator long() const throw(gmp_error);

Rational::operator int() const throw(gmp_error);

Rational::operator long() const throw(gmp_error);

Convert the value (the integral part of it) to the built-in type. If it is too large to fit into this type, raise an exception.

Integer::operator double() const;

Rational::operator double() const;

Find the next floating-point approximation.

std::swap(Integer&, Integer&);

std::swap(Rational&, Rational&);

Swap two values efficiently.

Integer& numerator(Rational&);

Integer& denominator(Rational&) throw(gmp_error);

const Integer& numerator(const Rational&);

const Integer& denominator(const Rational&) throw(gmp_error);

Give separate access to the numerator and denominator. Exception is raised by the attempt to assign 0 to the denominator.

We have intentionally not defined conversions from and to unsigned integral built-in types. It would have exploded the number of all possible function prototypes to an extent that would have been impossible to maintain. The unsigned numbers are almost never used for real arithmetic purposes, in opposite to the GMP numbers, so they shouldn't ever get in contact with each other. In the seldom case you really need it, you should cast the unsigned value to the signed type of the appropriate size.

Arithmetic operations

All arithmetic operators defined for the built-in integral types are available for every combination of Integer, Rational, int and long operands too. Surely, the assignment variants are also defined. The only exception makes up the residual operator %, which is not defined for rational numbers.

Unlike the complex classes in PTL, GMP wrappers don't implement reference counting nor lazy evaluation. Thus, all operators that seem to create an Integer or Rational object, are really doing this.

Integer& Integer::negate();

Rational& Rational::negate();

Efficient in-place negation.

int sign(const Integer&);

int sign(const Rational&);

Return an integer having the same sign as the given GMP number.

Integer std::abs(const Integer&);

Rational std::abs(const Rational&);

Return an absolute value.

Rational floor(const Rational&);

Rational ceil(const Rational&);

Find an integral number nearest to the given rational. The return type is declared Rational for the sake of congruence to the standard floor(double) and ceil(double) functions.

Integer div_exact(const Integer& a, const Integer& b);

Compute a/b efficiently, provided that b is a multiple of a. Otherwise result is undefined.

Integer::div_t std::div(const Integer& a, const Integer& b);

Compute a/b and a mod b simultaneously. The result type is built analogously to std::div_t:

struct Integer::div_t {
   Integer quot, rem;
};

Rational inv(const Rational& a);

Compute 1/a efficiently.

Integer pow(const Integer& a, unsigned long k);

Compute a^k.

Integer sqrt(const Integer& a);

Compute the integral part of square root.

Integer fac(unsigned long k);

Compute k! (factorial).

Integer binom(INT n, unsigned long k);

Compute the binomial coefficient n over k.

Integer gcd(INT a, INT b);

void gcd_ext(const Integer& a, const Integer& b, Integer& g, Integer& p, Integer& q);

Compute the greatest common divisor of two integral numbers. Extended algorithm also computes p and q such that g=a*p+b*q.

Integer lcm(INT a, INT b);

Compute the least common multiple of two integral numbers.

GMP also contains much more sofisticated number theoretic algorithms; we haven't needed them yet. They can be added to the interface by demand in a quite straightforward manner.

Comparisons

All comparison operators are defined for any combination of Integer, Rational, int and long operands.

int Integer::compare(const Integer& b) const;

int Rational::compare(const Rational& b) const;

Return -1 if *this is less than b, 0 if equal, 1 if greater.

Integer::operator bool() const;

Rational::operator bool() const;

bool Integer::operator! () const;

bool Rational::operator! () const;

Shortcut for an efficient comparison with 0 (useful in logical expressions).

bool abs_equal(const Number& a, const Number& b);

Efficient shortcut for abs(a)==abs(b)

namespace std_ext {
   template <> struct hash<Integer>;
   template <> struct hash<Rational>;
}

These specializations allow to use the GMP numbers as search keys in the containers std_ext::hash_set and std_ext::hash_map.

Input/output

explicit Integer(const char*) throw(gmp_error);

Integer::set(const char*) throw(gmp_error);

explicit Rational(const char*) throw(gmp_error);

Rational::set(const char*) throw(gmp_error);

Parse an ASCII representation and assign a converted value. The denominator must be separated by '/'; lacking denominator is assumed to be equal to 1. Numbers may start with 0 or 0x, which assumes octal or hexadecimal numeric base instead of the default decimal.

std::string Integer::to_string(int base=10);

std::string Rational::to_string(int base=10);

Creates an ASCII representation with a given numeric base.

std::istream& operator >> (std::istream&, Integer&);

std::istream& operator >> (std::istream&, Rational&);

Read a value from an input stream. The numeric base of the input ASCII representation is chosen according to the ios::basefield flags of the stream; if they are not set, the base is recognized automatically by the prefix: 0 (octal), 0x (hexadecimal), or default decimal. The input is consumed up to the first character not matching the chosen numeric base. If no characters were read at all, the ios::failbit of the stream is set.

std::ostream& operator << (std::ostream&, const Integer&);

std::ostream& operator << (std::ostream&, const Rational&);

Print the ASCII representation to the output stream. The numeric base and appearance are chosen according to the ios::basefield and ios::showbase flags of the stream; default base is 10. Rational numbers are printed in the form numerator/denominator; they are treated as one field with respect to the ios::witdh. By integral values the fractional line and denominator are suppressed.

Error handling

There are not so horribly many ways to play havoc with GMP numbers. They can't get overflowed as long as there is a piece of virtual memory available. The only remaining case is an attempt to divide thru 0. It is caught by GMP internally and mapped to a (platform-specific notion of) divide-by-zero signal. This is the same reaction one would expect from built-in types.

Unfortunately, GMP doesn't always test the denominators of Rationals; hence this test is performed in the wrapper class. Setting the denominator to zero will cause an exception of type

gmp_error : public std::domain_error

to be raised. Probably it would be better to emulate a zero division in this case too; every well-grounded opinion is welcome.

The second source of error conditions is parsing ASCII input. The handling of invalid data is context-dependent: input stream read operators report the error via the standard stream state flags, while the string conversion functions raise the gmp_error condition.