Step 1: Using the standard library

The basic data structures of the standard C++ library are available, once you have included the correct header files. Many of our constructions will be based on such containers. If you want for instance to manipulate lists of say integer, you will use the type:

  list<int>

which is a parametrised type list<T> in which we instantiate the type T by int. step1.cc Let us detail what is going on here:

1: As we use the list of the standard library, we include the header file list.h.
2: We define an abbreviation iterator for the iterators on the list of int.
6: We declare l as a list of int.
7: We append at the end of the list, the three elements 1,2,3.
9: We print it. We use an iterator to scan this list, from the beginning l.begin() to the end l.end(). See [STL96] for more information on these iterators and their correlation with generic programming.
10: We print a new line.

Once you have written this code in a file fich.cc, you can compile it with the usual compiler command

   g++ -o fich.ex fich.cc

and run the executable fich.ex or use the command alp:

   alp fich.cc

which will also produce the executable fich.ex. The result of this first execution should look like:

      1 2 3

For more information, on the basic data structures available in the stl, we refer to chapter X or to [STL96], [STL98].

Step 2: linear algebra

The main categories of objects for doing linear algebra are

VectStd: Standard dense vectors (X),
MatDense: Dense matrices (X),
MatStruct: Structured matrices (X),
MatSparse: Sparse matrices (X).

Before using these classes you will have to chose the container for the internal representation as illustrated now: If you want to construct standard vectors with double coefficients, you can use for instance the definition:

   typedef VectStd<array1d<double> > Vect;

The container array1d<C> (X) is a one-dimensional array of C, with a field specifying its size. But you can also used the containor vector<C> of the standard library or any other convenient container. See section X for more details. A vector can be constructed, for instance, from an array:

   double u[]={1,-2,1}; Vect V(3,u);

For the definition of dense matrices with double coefficients, you can use the definition:

   typedef MatDense<lapack<double> > Mat;

Here we could also have used the container array2d<C> (X) which extends array1d<C>, by specifying the number of rows and columns, but we prefer the lapack<C> cointainer (X) which will allow us to use the routines of the lapack library. The internal representation is also a one-dimensional array with a number of rows, a number of columns, a size bigger than the product of these two last numbers and a way to access to the element of index (i,j). There are many ways to build a matrix (see X). Here we parse a string:

   Mat A(3,3,"1 2 3 4 5 6 7 8 9 10");

It defines the following matrix:

      [ 1 2 3 ]
      [ 4 5 6 ]
      [ 7 8 9 ]

Check how are read the elements (in this construction it is by row). The usual operations on matrices and vectors are available:

   Vect W = A*V; cout <<W<<endl;

      [0,0,0]

   Mat B = A*A; cout <<B<<endl;

      [ 30 36 42 ]
      [ 66 81 96 ]
      [ 102 126 150 ]

   Mat C(-A); C += B + A*2; cout <<C<<endl;
 
      [ 31 38 45 ]
      [ 70 86 102 ]
      [ 109 134 159 ]

Views on vectors and matrices

These linear algebra classes also offers local views on their internal representation. A row (or a column) of a matrix, for instance, is a view on the corresponding elements of the matrix and share the properties of a vector. Indeed, it is a VectStd with a special container and can be manipulated as such. The types of the columns and the rows are accesible as

   MatDense<lapack<double> >::row_type
   MatDense<lapack<double> >::col_type

The characteristic of these views is that performing operations on the vectors will affect the object the view is looking at.

   Col(A,1)-= Col(A,0)*2;

transforms A into:

      [ 1  0 3 ]
      [ 4 -3 6 ]
      [ 7 -6 9 ]

It must be remembered that the memory space of the views Row(A,i), Col(A,i) is not duplicated but shared with the matrix A. It can be copied into a standard vector (each element if copied):

   W = Col(A,0);

Be carefull when you use these local views:

   Col(A,0)=Col(A,1);

means that the temporary object which is returned by Col(A,1) is assigned to the temporary object returned by Col(A,0). It does not mean that the coefficients viewed by Col(A,0) will be set to those viewed by Col(A,1). If this is what you want, here is the command to be used:

   Col(A,0)=Eval(Col(A,1));

Operations between rows (resp. columns) of different matrices are also possible:

   Mat D(2,2,"1 0 2 0");  Col(A,0)=Col(D,0)*10;

yields:

      [ 10  0 3 ]
      [ 20 -3 6 ]
      [ -6 -6 9 ]

Similarly subvectors also correspond to views on one-dimensional arrays. The operation

   W[Range(0,1)] = Eval(W[Range(1,2)]);

transforms the vector W=[1,4,7] into W=[4,4,7]. Or you can use operations such as

   Col(A,0)[Range(0,1)] = Eval(Col(A,1)[Range(1,2)]);

which transforms A into

      [ -3  0 3 ]
      [ -6 -3 6 ]
      [ -6 -6 9 ]

Here is the complete file for these examples, that you will compile with the command extttalp in order to get an executable file .ex: step2.cc

Step 3: univariate polynomials

Univariate polynomials are derived from the vectors. They have a multiplication, a degree and are printed using the variable x (X). Here is an example of definition:

   typedef UPolyDense<upar<Z<31> > > Pol;

We use the view of dense univariate polynomials UPolyDense, with the container upar which extends arra1d by storing the degree. The coefficients are of type Z<31> that is integers modulo 31 (X). In order to be able to use this class, we must first include the following files:

   #include "upoly.H"       // univariate polynomials
   #include "arithm/Zp.H"   // modular numbers

A polynomial can be constructed from an array of coefficient or parsed from a string (ending by a ;):

   Pol p("x^3-34*x+1;");

yields the polynomial

      (1)*x^3+(-3)*x+(23)

The ``natural'' operations of polynomials are available (see X):

   Pol q = p*p;

yields the polynomial

      (1)*x^6+(-6)*x^4+(15)*x^3+(9)*x^2+(-14)*x+(2)

Here all the operations on the coefficients are performed modulo 31. These polynomials also have an euclidean division:

   p+= Pol("x+1;");  Pol r= q%p;

yields the remainder r

   (1)*x^2+(-29)*x+(1)

which is also (x+1)² mod 31. Notice that the operator % will yield the remainder of the Euclidean division if the leading term is invertible (here it is the case for 31 is prime). The quotient is given by the operator /. Univariate polynomials can also be evaluated:

   cout <<q(Z<31>(1))<<endl;

yields 9, which is of type Z<31>. The output of the evaluation is of the same type as the input type, if it is compatible with the coefficient type. For instance, q(1) will not be valid here for the evaluation tries to output an int from arithmetic operations between int and Z<31>¹ .

Next is another example where we use big integers, defined from the file bignum.h (X):

  typedef UPolyDense<upar<Scl<MPZ> > >   ZPol;
  ZPol a("x^2-1243432232*x+2334345435434;");
  ZPol b=a^3;    cout <<b<<endl;

     (1)*x^6+(-3730296696)*x^5
     +(4638378149765811774)*x^4
     +(-1922517478209553074031443296)*x^3
     +(10827576861722625150704513999916)*x^2
     +(-20327015669035171107600087000151776)*x
     +(12720241876172641933421003452296326504)

Here is the complete file for this section: step3.cc

Step 4: multivariate polynomials

Multivariate monomials are defined by a type C for the coefficients and a type E for the exponents (X):

   Monom<C,E>

Multivariate polynomials are represented as sorted list of monomials, defined by a container R and a trait class type O for the order on the monomials (X):

   MPoly<R,O>

Let us detail their use on the following example (file step4.cc): step4.cc

1: Inclusion of the necessary header files, for the manipulation of multivariate polynomials.
3: We define an abbreviation for the data-structure, corresponding to monomials. These are monomials with double coefficients and dynamicexp<'x'> for the exponent type (see X).
4: We define an abbreviation for the data-structure, corresponding to the polynomials. They are represented by a list of monomials, sorted by Degree Lexicographic ordering (see X).
6:We define the main procedure.
8: We build the polynomial p, by parsing the string "x[1]^2*x[2]+x[1]+1". Other techniques for constructing p are detailed in section X.
10: We multiply p by q.
11: We output the new value of the variable p.
12: We scan the monomials of p from the beginning to the end.
13: And we output the coefficients of the monomials.

Once you have typed this file, you compile it with the command alp, run it and obtain the following result:

   x1^2*x2+(2)*x1+(3)
   x1^4*x2+x1^2*x2*x3+(2)*x1^3+(3)*x1^2+(2)*x1*x3+(3)*x3
   1 1 2 3 2 3

Troubleshooting

Things which are not allowed for the moment:

Avoid
```
 p = p*p; p*=p; b = A*b; ..
```
Instead, use
```
 p = Pol(p*p); p *= Pol(p); b = Vect(A*b); ..
```
which make a copy of the right hand side, in order to avoid problems when assigning a value to the left hand side, which could depend on this left hand side.
...

1: In a next version, we plan to implement traits classes which will simplify these problems.