C++ For Mathematicians (2006) [eng]

The GMP package provides four important C++ classes:

•mpz_class for arbitrary precision integers,

•mpq_class for exact rational numbers,

•mpf_class for floating point numbers, and

•gmp_randclass for generating large random numbers.

To use these classes, your program needs to include a header file:

#include <gmpxx.h>

Variables are declared in the usual manner. For example, to declare x to be an arbitrary precision integer, use this:

mpz_class x;

This initializes x with the value 0. To give it a large value, the following does not work.

The problem is that C++ is not able to deal with that large value as one of its basic types. Instead, you can type this:

x = "3098472938750987439857234543534232626245985"; // this works

The mpz_class type can convert character arrays (containing digits). The mpz_class can be initialized using other bases. For example,

mpz_class x("12321423112312321312314001200001213", 5);

initializes x with a base-5 value.

The usual arithmetic and comparison operators work just as expected. Multiplication of extremely large integers is accomplished using sophisticated efficient algorithms.

The rational type, mpq_class, is constructed either from two integer arguments

mpq_class a(n_1,n_2);

or else a character array, such as this:

mpq_class b("53490875234097/1134381");

Rationals can also be constructed from double values. Indeed, the GMP types can be converted to and constructed from nearly any numeric type.

Rational mpq_class objects have a canonicalize() method. The purpose of this method is to clear any common factors between numerator and denominator. It’s a good idea to invoke this method after creating an mpq_class value.

The C++ features of the GMP package are a supplement to the C base that forms the bulk of GMP. Some GMP procedures require some fancy footwork to be used in C++. For example, to find the greatest common divisor of two mpz_class integers one uses the procedure mpz_gcd. This procedure takes three arguments of type mpz_t—not mpz_class. The mpz_class objects contain an mpz_t type value internally and to access that internal value one uses the get_mpz_t() method. Here’s how this all works.

First, we set up our variables:

mpz_class a = "47825100"; mpz_class b = "55431225";

mpz_class d; // place to hold the answer

Then we call the mpz_gcd method like this:

mpz_gcd(d.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());

Now d holds the gcd of a and b.

An object of type gmp_randclass is a random number generator. This object offers a choice of pseudo random number generator algorithms and allows the user to set the seed. For example, to create a new random number generator with a default algorithm and seeded from the system clock, we write this:

gmp_randclass X(gmp_randinit_default); X.seed(time(0));

To extract a random value from the generator, we use the get_z_bits method and specify the size (in bits) of the result. For example, X.get_z_bits(100) returns a 100-bit random integer.

Here is a program that illustrates these ideas.

Program 13.1: A program to illustrate the use of the GMP package.

1 #include <gmpxx.h>

2#include <iostream>

3

4using namespace std;

5

6int main() {

7mpz_class a, b;

8

9 a = "54098745908347598037452";

10 b = "44523409864";

11

12cout << "a = " << a << endl;

13cout << "b = " << b << endl;

14cout << "a*b = " << a*b << endl;

15cout << "a/b = " << a/b << endl;

16cout << "a%b = " << a%b << endl;

17cout << "a+b = " << a+b << endl;

18cout << "a-b = " << a-b << endl;

19cout << "b-a = " << b-a << endl;

On my computer, the gmpxx.h header file is in the directory /sw/include (specified by the -I option), the libraries are located in /sw/lib (specified by the -L option), and the gmp and gmpxx libraries are named explicitly (by the pair of -l options).

13.2Linear algebra

13.2.1 Two-dimensional arrays in C++

One-dimensional arrays in C++ are mildly difficult to use. One needs either to give the explicit size of the array in advance (e.g., long vals[10];) or else to declare the array via a pointer (long *vals;), allocate space (vals = new long[n];), and then release the memory when finished (delete[] vals;). We need to remember that the first element of an array is indexed as vals[0] and there is no protection against accessing data beyond the bounds of the array. There is no convenient way to extend the size of the array. The vector class template solves many of these problems.

It is possible to use multidimensional arrays in C++, but the difficulties are compounded. We describe such arrays here briefly in order to convince you that you do not want to use these tricky structures and, instead, will avail yourself of some of the options we discuss below.

To declare a two-dimensional array of a given size, we use a statement such as this:

long vals[4][10];

This establishes a variable named vals as a 4 × 10-array of long integer values. Please note that the following statement is incorrect: long vals[4,10];. Unfortunately, it is valid C++ (don’t bother trying to figure out what it does) so the compiler won’t catch this error.

To access an element of this array, we use the syntax vals[i][j] where i is between 0 and 3, and j is between 0 and 9. The expression vals[i,j] is incorrect.

It is possible to declare three- (or higher-) dimensional arrays. For example,

long vals[4][10][19];

declares vals to be a 4 ×10 ×19-array of long values. The i, j,k-entry of this array is given by vals[i][j][k] where i, j,k have the expected constraints.

Passing multidimensional arrays to procedures is possible, but the syntax is tricky. Procedures that receive such arrays need to specify the size of the array they expect to receive.

Fortunately, there are better alternatives. The two that we present have the added advantage of providing linear algebraic capabilities (matrix multiplication, eigenvalue calculation, etc.).

13.2.2 The TNT and JAMA packages

The U.S. government’s National Institute of Standards and Technology (NIST) has created a pair of packages for working with one-, two-, and three-dimensional arrays and to perform standard linear algebra functions thereon. The packages are the Template Numerical Toolkit (TNT) and the C++ Java Matrix Library (JAMA).2

These are available for free download from the following Web site,

http://math.nist.gov/tnt/index.html

You need to download the TNT and JAMA packages separately. Be sure to download the documentation as well; this consists of a collection of Web pages that describe the various classes and procedures. (We describe TNT version 1.2.4 and JAMA version 1.2.2.)

Installation is easy because these packages are simply collections of .h header files. They can be copied to any convenient location on your computer (including, e.g., the directory containing your program that uses these files).

The most important classes in the TNT package are Array1D (for vectors) and Array2D (for matrices3). To use these, we need the directive #include "tnt.h" and the statement using namespace TNT; at the beginning of our program. The classes and procedures in the TNT class are in their own namespace. If we did not use this statement, then we would need to prepend TNT:: to the names of TNT classes and procedures (e.g., TNT::Array2D in lieu of Array2D).

To declare an array of elements of type, say, double, we use statements like these:

Array2D<double> M(5,10);

Array1D<double> v(9);

Array2D<double> X;

2This library was developed first for the Java programming language and subsequently ported to C++.

3You may discover that the TNT package includes classes called Matrix and Vector. Do not use these. They are deprecated classes; this means they are obsolete and their inclusion in the package is only to support people who used older versions of TNT.

Here, M is a 5 × 10-array (matrix) of double values, v is a length-9 array (vector), and X is an as-yet unsized matrix.

To access elements of an array we use the usual C++ conventions. For the variables M and v above, the syntax is as follows,

•M[i][j] where 0 ≤ i < 5 and 0 ≤ j < 10, and

•v[k] where 0 ≤ k < 9.

To learn the size of an array, the methods dim1() and dim2() give the number of rows and columns (if more than one-dimensional) in the array.

Read this carefully: The assignment operator for TNT array classes (e.g., B=A;) works in a manner that is different from other classes we have encountered thus far. After this assignment not only are A and B arrays of the same size and not only do they hold the same values, they now refer to the exact same data. In other words, modifications to one of A or B results in change in the other. (We say that A and B provide different views to the same underlying data.) If we wish to make an independent copy of A in B, we need to use the copy method like this:

B = A.copy();

Here is a program that illustrates the difference between the statements B=A; and

B=A.copy();.

Program 13.2: Assignment versus copying in the TNT package.

1 #include <iostream>

2#include "tnt.h"

3

4 using namespace std;

5using namespace TNT;

6

7int main() {

8Array2D<double> A(3,3), B, C;

9

10for (int i=0; i<3; i++)

11for (int j=0; j<3; j++)

12A[i][j] = (i+1) + 0.1*(j+1);

19B[1][1] = 888;

20A[0][2] = 999;

21

22cout << "Now A is this " << endl << A << endl;

23cout << "and B is this " << endl << B << endl;

24cout << "and C is this " << endl << C << endl;

25

26return 0;

27}