Snippets of code

Simple arithmetic with polynomials

#include "my_gmp.H"
#include "algebras.H"
#include <iostream>
using namespace std;


int main(){
	AffineAlgebra<ZZ> A;
	Polynomial<ZZ> I, x, y, P, Q;

	x= A.new_variable("X");
	y= A.new_variable("Y");
	I= A.one();
	
	P= x*x + y*x + y*5;
	Q= I + x*(-1) + y*y*y;

	cout << "P= " << P << endl
		<< "Q= " << Q << endl
		<< "and P*Q=" << P*Q << endl; 	


	return 0;
}

P= X^2 + 5Y + XY
Q= I + -1X + Y^3
and P*Q=X^2 + -1X^3 + 5Y + -4XY + -1X^2Y + X^2Y^3 + 5Y^4 + XY^4

• I use the GMP "bignum" library. Here the my_gmp.H file which gets included is simply here to give new names to the GMP objects and functions: notably, ZZ for the integers, and QQ for the rationals (there is also the gcd function for well, the gcd). These replace the original mpz_int, etc. Later you will also find the type F_2 for the field with two elements.
• The algebras.H contains everything about algebras, polynomials, modules and so on.
• We first create an affine (=finitely generated) algebra over the integers ZZ. Originally it is just equal to ZZ, and after we add two variables, it is a polynomial algebra ZZ[X,Y]. Notice that the function which creates the new variables also has a value, of type Polynomial< ZZ >, which we keep for later computations. Also, we ask A to give us its unit with the one() method.
• The variable x could have been given any name, of course; it refers to the element X of the algebra A. Similary for y and Y. Their use is pretty intuitive. Notice that we write y*5, while 5*y would create an error. Also, we have written Q= I + ..., not just Q= 1 + ..., which would also create an error. The polynomial X - 2 can be created with x + I*(-2).

Quotient algebras

#include "my_gmp.H"
#include "algebras.H"
#include <iostream>
#include "ftwo.H"
using namespace std;


int main(){
	AffineAlgebra<F_2> A;
	Polynomial<F_2> I, u, x;
	//let us create the field with 4 elements:
	u= A.new_variable("u");
	I= A.one();
	A.add_relation(I + u + u*u);
	
	x= I+u;
	
	cout << "any element to the third power is 1..." << endl
		<< "and indeed (1+u)^3= " << x*x*x << endl
		<< "which is also " << A.reduced_form(x*x*x) << endl;

	return 0;
}

any element to the third power is 1...
and indeed (1+u)^3= I + u + u^2 + u^3
which is also I

• Here the method add_relation creates a quotient algebra, namely A is now divided by the ideal generated by the relation given. You can add any number of relations.
• The method reduced_form performs a division algorithm -- when there is only one relation, as in the example, it is just Euclidean division. However, when more relations are given, this will only work well after a Grobner basis has been computed: see the next example.

Grobner basis

#include "my_gmp.H"
#include "algebras.H"
#include <iostream>
#include "ftwo.H"
using namespace std;


int main(){
	AffineAlgebra<F_2> A;
	Polynomial<F_2> I, y, x;
	x= A.new_variable("x");
	y= A.new_variable("y");
	I= A.one();
	A.add_relation(x*x + y*y +x*y);
	A.add_relation(x*x*y + x*y*y);
	
	cout << A << endl;
	cout << "y^4=" << A.reduced_form(y*y*y*y) << endl;
	
	A.find_grobner_basis();
	
	cout << endl << A << endl;
	cout << "y^4=" << A.reduced_form(y*y*y*y) << endl;
	
	return 0;
}

generators:
x y 
relations:
x^2 + xy + y^2
x^2y + xy^2

y^4=x^3y

generators:
x y 
relations:
x^2 + xy + y^2
x^3

y^4=0

• We create an algebra with two generators and two relations. Notice how, when we ask for the reduced form of y^4, we get a wrong (= not reduced) answer (though the equality y^4 = x^3y does hold, of course!).
• After find_grobner_basis() has been called, the relations have changed in A, they now form a Grobner basis, and the reduced form obtained is correct.
• The order on the variables (needed when you talk about Grobner basis) is chronological: here x < y for example. The order on power products is LEX by default. You can change this behaviour by adding the line

  PowProd::use_order(DEGLEX);

  at the beginning of the program. DEGREVLEX and a couple of others are available.
  Warning: one should pick an order and stick to it. Changing the order will screw up the objects (polynomials and algebras) already created. (Just in case anyone were trying to experiment with the code, which I doubt).

Do algebras have elements?

#include "my_gmp.H"
#include "algebras.H"
#include <iostream>
#include "ftwo.H"
using namespace std;


int main(){
	AffineAlgebra<ZZ> A;
	Polynomial<ZZ> I, y, x;
	x= A.new_variable("x");
	y= A.new_variable("y");
	
	AffineAlgebra<ZZ> B;
	Polynomial<ZZ> J, u, v;
	u= B.new_variable("u");
	v= B.new_variable("v");
	
	Polynomial<ZZ> P;
	P= x*x + y*4 + x*y*2;
	
	cout << "P= " << P << endl;
	P.alphabet= &B;
	cout << "and now P= " << P << endl;
	
	return 0;
}

P= x^2 + 4y + 2xy
and now P= u^2 + 4v + 2uv

• As you can see on this example, a polynomial has an attribute "alphabet", which is a pointer to an algebra (most of the time -- it could also point to an abstract class called Alphabet, from which algebras are derived). It is used only for display purposes. You can see that, by changing the alphabet, the polynomial prints differently.
• On this example, the alphabet of P is originally A, because it is so for x and y, which were created using A.
• Most of the time, one should avoid switching the alphabet like this: it is dangerous, for a crash can occur in the situation where B has less variables than A.

Copying algebras

	B= A;
	B.fix_alphabets();

• The command B= A; copies everything in A and puts it in B. Which means that the relations in A, if any, are copied as such, and they are polynomials whose alphabets points at A. Now, A could disappear from the memory while B lingers, of course, so it is desirable to have the relations in B have the "right" alphabet, namely B. This is what fix_alphabets() does.
• Of course, I should have overloaded the operator = so one does not have to call fix_alphabets(). Since there are quite a few types of algebras, I was too lazy to do so (the original operator = gets most of the job right!). I guess I should at least have written a copy() method. Oh well.