8.6 Numerics

The C++ library has several headers that support numerical programming. The most basic header is <cmath>, which declares transcendental and other mathematical functions. This C header is expanded in C++ to declare overloaded versions of every function. For example, whereas C declares only exp(double), C++ also declares exp(float) and exp(long double).

The <complex> header declares a class template for complex numbers, with the specializations you would probably expect for float, double, and long double. Transcendental and I/O functions are also declared for complex numbers.

The <numeric> header declares a few algorithms (that use standard iterators) for numerical sequences.

The most interesting numerical functions are in <valarray>. A valarray is like an ordinary numerical array, but the compiler is free to make some simplifying assumptions to improve optimization. A valarray is not a container, so it cannot be used with standard algorithms or the standard numeric algorithms.

You can use the complex template as an example of how to define custom numeric types. For example, suppose you want to define a type to represent rational numbers (fractions). To use rational objects in a valarray, you must define the class so it behaves the same as ordinary values, such as ints. In other words, a custom numeric type should have the following:

• A public default constructor (e.g., rational( ))

• A public copy constructor (e.g., rational(const rational&))

• A public destructor

• A public assignment operator (e.g., rational& operator=(const rational&))

• Reasonable arithmetic and comparison operators

When you implement a class, make sure that the copy constructor and assignment operator have similar, reasonable results.

When overloading arithmetic and comparison operators, think about which operators are meaningful (e.g., most arithmetic types should have addition, subtraction, multiplication, and division, but only integral types will probably have remainder, bitwise, and shift operators). You should provide overloaded operators that accept the built-in types as operands. For example, the complex template defines the following functions for operator+:

template<typename T> complex<T> operator+(const complex<T>& z);
template<typename T> complex<T> operator+(const complex<T>& x, const
                                          complex<T>& y);
template<typename T> complex<T> operator+(const complex<T>& x, const T& y);
template<typename T> complex<T> operator+(const T& x, const complex<T>& y);

Example 8-5 shows excerpts from a simple rational class template.

Example 8-5. The rational class template for rational numbers

template<typename T>
class rational
{
public:
  typedef T value_type;

  rational(  )                 : num_(0),   den_(1) {}
  rational(value_type num)   : num_(num), den_(1) {}
  rational(value_type num, value_type den)
    : num_(num), den_(den) { reduce(  ); }
  rational(const rational& r): num_(r.num_), den_(r.den_) {}
  template<typename U>
  rational(const rational<U>& r)
    : num_(r.num_), den_(r.den_) { reduce(  ); }

  rational& operator=(const rational& r)
    { num_ = r.num_; den_ = r.den_; return *this; }
  template<typename U>
  rational& operator=(const rational<U>& r)
    { assign(r.numerator(  ), r.denominator(  )); return *this; }

  void assign(value_type n, value_type d)
    { num_ = n; den_ = d; reduce(  ); }

  value_type numerator(  )   const { return num_; }
  value_type denominator(  ) const { return den_; }

private:
  void reduce(  );
  value_type num_;
  value_type den_;
};

// Reduce the numerator and denominator by the gcd. Make sure that the
// denominator is nonnegative.
template<typename T>
void rational<T>::reduce(  )
{
  if (den_ < 0) {
    den_ = -den_;
    num_ = -num_;
  }
  T d = gcd(num_, den_);
  num_ /= d;
  den_ /= d;
}

// Greatest common divisor using Euclid's algorithm
template<typename T>
T gcd(T n, T d)
{
  n = abs(n);
  while (d != 0) {
    T t = n % d;
    n = d;
    d = t;
  }
  return n;
}

// Multiplication assignment operator. Often implemented as a member function,
// but there is no need to do so.
template<typename T, typename U>
rational<T>& operator*=(rational<T>& dst,
                        const rational<U>& src)
{
  dst.assign(dst.numerator(  ) * src.numerator(  ),
             dst.denominator(  ) * src.denominator(  ));
  return dst;
}

// Multiply two rational numbers.
template<typename T>
rational<T> operator*(const rational<T>& a,
                      const rational<T>& b)
{
  rational<T> result(a);
  result *= b;
  return result;
}

// Multiply rational by an integral value.
template<typename T>
rational<T> operator*(const T& a, const rational<T>& b)
{
  return rational<T>(a * b.numerator(  ), b.denominator(  ));
}

template<typename T>
rational<T> operator*(const rational<T>& a, const T& b)
{
  return rational<T>(b * a.numerator(  ), a.denominator(  ));
}
// Other arithmetic operators are similar.
// Comparison. All other comparisons can be implemented in terms of operator==
// and operator<.
template<typename T>
bool operator==(const rational<T>& a, const rational<T>& b)
{
  // Precondition. Both operands are reduced.
  return a.numerator(  ) == b.numerator(  ) &&
         a.denominator(  ) == b.denominator(  );
}

template<typename T>
bool operator<(const rational<T>& a, const rational<T>& b)
{
  return a.numerator(  ) * b.denominator(  ) <
         b.numerator(  ) * a.denominator(  );
}

Many numerical programmers find the C++ standard library to be lacking. However, the Blitz++ project is a popular, high-performance numerical library. Boost also has some numerical headers, such as a full-featured rational class template. See Appendix B for information about these and other C++ libraries.