In this section we will show how one can add a new domain to GAP. All domains are implemented in the library in this way. We will use the ring of Gaussian integers as our example.
Note that everything defined here is already in the library file
LIBNAME/"gaussian.g", so there is no need for you to type it in. You
may however like to make a copy of this file and modify it.
The elements of this domain are already available, because Gaussian
integers are just a special case of cyclotomic numbers. As is described
in chapter Cyclotomics E(4) is GAP's name for the complex root of
-1. So all Gaussian integers can be represented as a + b*E(4),
where a and b are ordinary integers.
As was already mentioned each domain is represented by a record. So we
create a record to represent the Gaussian integers, which we call
GaussianIntegers.
gap> GaussianIntegers := rec();;
The first components that this record must have are those that identify this record as a record denoting a ring domain. Those components are called the category components.
gap> GaussianIntegers.isDomain := true;;
gap> GaussianIntegers.isRing := true;;
The next components are those that uniquely identify this ring. For
rings this must be generators, zero, and one. Those components are
called the identification components of the domain record. We also
assign a name component. This name will be printed when the domain is
printed.
gap> GaussianIntegers.generators := [ 1, E(4) ];;
gap> GaussianIntegers.zero := 0;;
gap> GaussianIntegers.one := 1;;
gap> GaussianIntegers.name := "GaussianIntegers";;
Next we enter some components that represent knowledge that we have about
this domain. Those components are called the knowledge components. In
our example we know that the Gaussian integers form a infinite,
commutative, integral, Euclidean ring, which has an unique factorization
property, with the four units 1, -1, E(4), and -E(4).
gap> GaussianIntegers.size := "infinity";;
gap> GaussianIntegers.isFinite := false;;
gap> GaussianIntegers.isCommutativeRing := true;;
gap> GaussianIntegers.isIntegralRing := true;;
gap> GaussianIntegers.isUniqueFactorizationRing := true;;
gap> GaussianIntegers.isEuclideanRing := true;;
gap> GaussianIntegers.units := [1,-1,E(4),-E(4)];;
This was the easy part of this example. Now we have to add an
operations record to the domain record. This operations record
(GaussianIntegers.operations) shall contain functions that implement
all the functions mentioned in chapter Rings, e.g., DefaultRing,
IsCommutativeRing, Gcd, or QuotientRemainder.
Luckily we do not have to implement all this functions. The first class
of functions that we need not implement are those that can simply get the
result from the knowledge components. E.g., IsCommutativeRing looks
for the knowledge component isCommutativeRing, finds it and returns
this value. So GaussianIntegers.operations.IsCommutativeRing is never
called.
gap> IsCommutativeRing( GaussianIntegers );
true
gap> Units( GaussianIntegers );
[ 1, -1, E(4), -E(4) ]
The second class of functions that we need not implement are those for which there is a general algorithm that can be applied for all rings. For example once we can do a division with remainder (which we will have to implement) we can use the general Euclidean algorithm to compute the greatest common divisor of elements.
So the question is, how do we get those general functions into our
operations record. This is very simple, we just initialize the
operations record as a copy of the record RingOps, which contains all
those general functions. We say that GaussianIntegers.operations
inherits the general functions from RingOps.
gap> GaussianIntegersOps := OperationsRecord(
> "GaussianIntegersOps", RingOps );;
gap> GaussianIntegers.operations := GaussianIntegersOps;;
So now we have to add those functions whose result can not (easily) be
derived from the knowledge components and that we can not inherit from
RingOps.
The first such function is the membership test. This function must test
whether an object is an element of the domain GaussianIntegers.
IsCycInt(x) tests whether x is a cyclotomic integer and
NofCyc(x) returns the smallest n such that the cyclotomic x can
be written as a linear combination of powers of the primitive n-th root
of unity E(n). If NofCyc(x) returns 1, x is an ordinary
rational number.
gap> GaussianIntegersOps.\in := function ( x, GaussInt )
> return IsCycInt( x ) and (NofCyc( x ) = 1 or NofCyc( x ) = 4);
> end;;
Note that the second argument GaussInt is not used in the function.
Whenever this function is called, the second argument must be
GaussianIntegers, because GaussianIntegers is the only domain that
has this particular function in its operations record. This also happens
for most other functions that we will write. This argument can not be
dropped though, because there are other domains that share a common in
function, for example all permutation groups have the same in function.
If the operator in would not pass the second argument, this function
could not know for which permutation group it should perform the
membership test.
So now we can test whether a certain object is a Gaussian integer or not.
gap> E(4) in GaussianIntegers;
true
gap> 1/2 in GaussianIntegers;
false
gap> GaussianIntegers in GaussianIntegers;
false
Another function that is just as easy is the function Random that
should return a random Gaussian integer.
gap> GaussianIntegersOps.Random := function ( GaussInt )
> return Random( Integers ) + Random( Integers ) * E( 4 );
> end;;
Note that actually a Random function was inherited from RingOps. But
this function can not be used. It tries to construct the sorted list of
all elements of the domain and then picks a random element from that
list. Therefor this function is only applicable for finite domains, and
can not be used for GaussianIntegers. So we overlay this default
function by simply putting another function in the operations record.
Now we can already test whether a Gaussian integer is a unit or not.
This is because the default function inherited from RingOps tests
whether the knowledge component units is present, and it returns true
if the element is in that list and false otherwise.
gap> IsUnit( GaussianIntegers, E(4) );
true
gap> IsUnit( GaussianIntegers, 1 + E(4) );
false
Now we finally come to more interesting stuff. The function Quotient
should return the quotient of its two arguments x and y. If the
quotient does not exist in the ring (i.e., if it is a proper Gaussian
rational), it must return false. (Without this last requirement we
could do without the Quotient function and always simply use the /
operator.)
gap> GaussianIntegersOps.Quotient := function ( GaussInt, x, y )
> local q;
> q := x / y;
> if not IsCycInt( q ) then
> q := false;
> fi;
> return q;
> end;;
The next function is used to test if two elements are associate in the
ring of Gaussian integers. In fact we need not implement this because
the function that we inherit from RingOps will do fine. The following
function is a little bit faster though that the inherited one.
gap> GaussianIntegersOps.IsAssociated := function ( GaussInt, x, y )
> return x = y or x = -y or x = E(4)*y or x = -E(4)*y;
> end;;
We must however implement the function StandardAssociate. It should
return an associate that is in some way standard. That means, whenever
we apply StandardAssociate to two associated elements we must obtain
the same value. For Gaussian integers we return that associate that lies
in the first quadrant of the complex plane. That is, the result is that
associated element that has positive real part and nonnegative imaginary
part. 0 is its own standard associate of course. Note that this is a
generalization of the absolute value function, which is
StandardAssociate for the integers. The reason that we must implement
StandardAssociate is of course that there is no general way to compute
a standard associate for an arbitrary ring, there is not even a standard
way to define this!
gap> GaussianIntegersOps.StandardAssociate := function ( GaussInt, x )
> if IsRat(x) and 0 <= x then
> return x;
> elif IsRat(x) then
> return -x;
> elif 0 < COEFFSCYC(x)[1] and 0 <= COEFFSCYC(x)[2] then
> return x;
> elif COEFFSCYC(x)[1] <= 0 and 0 < COEFFSCYC(x)[2] then
> return - E(4) * x;
> elif COEFFSCYC(x)[1] < 0 and COEFFSCYC(x)[2] <= 0 then
> return - x;
> else
> return E(4) * x;
> fi;
> end;;
Note that COEFFSCYC is an internal function that returns the
coefficients of a Gaussian integer (actually of an arbitrary cyclotomic)
as a list.
Now we have implemented all functions that are necessary to view the Gaussian integers plainly as a ring. Of course there is not much we can do with such a plain ring, we can compute with its elements and can do a few things that are related to the group of units.
gap> Quotient( GaussianIntegers, 2, 1+E(4) );
1-E(4)
gap> Quotient( GaussianIntegers, 3, 1+E(4) );
false
gap> IsAssociated( GaussianIntegers, 1+E(4), 1-E(4) );
true
gap> StandardAssociate( GaussianIntegers, 3 - E(4) );
1+3*E(4)
The remaining functions are related to the fact that the Gaussian integers are an Euclidean ring (and thus also a unique factorization ring).
The first such function is EuclideanDegree. In our example the
Euclidean degree of a Gaussian integer is of course simply its norm.
Just as with StandardAssociate we must implement this function because
there is no general way to compute the Euclidean degree for an arbitrary
Euclidean ring. The function itself is again very simple. The Euclidean
degree of a Gaussian integer x is the product of x with its complex
conjugate, which is denoted in GAP by GaloisCyc( x, -1 ).
gap> GaussianIntegersOps.EuclideanDegree := function ( GaussInt, x )
> return x * GaloisCyc( x, -1 );
> end;;
Once we have defined the Euclidean degree we want to implement the
QuotientRemainder function that gives us the Euclidean quotient and
remainder of a division.
gap> GaussianIntegersOps.QuotientRemainder := function ( GaussInt, x, y )
> return [ RoundCyc( x/y ), x - RoundCyc( x/y ) * y ];
> end;;
Note that in the definition of QuotientRemainder we must use the
function RoundCyc, which views the Gaussian rational x/y as a
point in the complex plane and returns the point of the lattice spanned
by 1 and E(4) closest to the point x/y. If we would truncate
towards the origin instead (this is done by the function IntCyc) we
could not guarantee that the result of EuclideanRemainder always has
Euclidean degree less than the Euclidean degree of y as the following
example shows.
gap> x := 2 - E(4);; EuclideanDegree( GaussianIntegers, x );
5
gap> y := 2 + E(4);; EuclideanDegree( GaussianIntegers, y );
5
gap> q := x / y; q := IntCyc( q );
3/5-4/5*E(4)
0
gap> EuclideanDegree( GaussianIntegers, x - q * y );
5
Now that we have implemented the QuotientRemainder function we can
compute greatest common divisors in the ring of Gaussian integers. This
is because we have inherited from RingOps the general function Gcd
that computes the greatest common divisor using Euclid's algorithm,
which only uses QuotientRemainder (and StandardAssociate to return
the result in a normal form). Of course we can now also compute least
common multiples, because that only uses Gcd.
gap> Gcd( GaussianIntegers, 2, 5 - E(4) );
1+E(4)
gap> Lcm( GaussianIntegers, 2, 5 - E(4) );
6+4*E(4)
Since the Gaussian integers are a Euclidean ring they are also a unique factorization ring. The next two functions implement the necessary operations. The first is the test for primality. A rational integer is a prime in the ring of Gaussian integers if and only if it is congruent to 3 modulo 4 (the other rational integer primes split into two irreducibles), and a Gaussian integer that is not a rational integer is a prime if its norm is a rational integer prime.
gap> GaussianIntegersOps.IsPrime := function ( GaussInt, x )
> if IsInt( x ) then
> return x mod 4 = 3 and IsPrimeInt( x );
> else
> return IsPrimeInt( x * GaloisCyc( x, -1 ) );
> fi;
> end;;
The factorization is based on the same observation. We compute the
Euclidean degree of the number that we want to factor, and factor this
rational integer. Then for every rational integer prime that is
congruent to 3 modulo 4 we get one factor, and we split the other
rational integer primes using the function TwoSquares and test which
irreducible divides.
gap> GaussianIntegersOps.Factors := function ( GaussInt, x )
> local facs, # factors (result)
> prm, # prime factors of the norm
> tsq; # representation of prm as x^2 + y^2
>
> # handle trivial cases
> if x in [ 0, 1, -1, E(4), -E(4) ] then
> return [ x ];
> fi;
>
> # loop over all factors of the norm of x
> facs := [];
> for prm in Set( FactorsInt( EuclideanDegree( x ) ) ) do
>
> # p = 2 and primes p = 1 mod 4 split according to p = x^2+y^2
> if prm = 2 or prm mod 4 = 1 then
> tsq := TwoSquares( prm );
> while IsCycInt( x / (tsq[1]+tsq[2]*E(4)) ) do
> Add( facs, (tsq[1]+tsq[2]*E(4)) );
> x := x / (tsq[1]+tsq[2]*E(4));
> od;
> while IsCycInt( x / (tsq[2]+tsq[1]*E(4)) ) do
> Add( facs, (tsq[2]+tsq[1]*E(4)) );
> x := x / (tsq[2]+tsq[1]*E(4));
> od;
>
> # primes p = 3 mod 4 stay prime
> else
> while IsCycInt( x / prm ) do
> Add( facs, prm );
> x := x / prm;
> od;
> fi;
>
> od;
>
> # the first factor takes the unit
> facs[1] := x * facs[1];
>
> # return the result
> return facs;
> end;;
So now we can factorize numbers in the ring of Gaussian integers.
gap> Factors( GaussianIntegers, 10 );
[ -1-E(4), 1+E(4), 1+2*E(4), 2+E(4) ]
gap> Factors( GaussianIntegers, 103 );
[ 103 ]
Now we have written all the functions for the operations record that
implement the operations. We would like one more thing however. Namely
that we can simply write Gcd( 2, 5 - E(4) ) without having to specify
GaussianIntegers as first argument. Gcd and the other functions
should be clever enough to find out that the arguments are Gaussian
integers and call GaussianIntegers.operations.Gcd automatically.
To do this we must first understand what happens when Gcd is called
without a ring as first argument. For an example suppose that we have
called Gcd( 66, 123 ) (and want to compute the gcd over the integers).
First Gcd calls DefaultRing( [ 66, 123 ] ), to obtain a ring that
contains 66 and 123. DefaultRing then calls Domain( [ 66, 123 ] ) to
obtain a domain, which need not be a ring, that contains 66 and 123.
Domain is the only function in the whole GAP library that knows
about the various types of elements. So it looks at its argument and
decides to return the domain Integers (which is in fact already a ring,
but it could in principle also return Rationals). DefaultRing now
calls Integers.operations.DefaultRing( [ 66, 123 ] ) and expects a ring
in which the requested gcd computation can be performed.
Integers.operations.DefaultRing( [ 66, 123 ] ) also returns Integers.
So DefaultRing returns Integers to Gcd and Gcd finally calls
Integers.operations.Gcd( Integers, 66, 123 ).
So the first thing we must do is to tell Domain about Gaussian
integers. We do this by extending Domain with the two lines
elif ForAll( elms, IsGaussInt ) then
return GaussianIntegers;
so that it now looks as follows.
gap> Domain := function ( elms )
> local elm;
> if ForAll( elms, IsInt ) then
> return Integers;
> elif ForAll( elms, IsRat ) then
> return Rationals;
> elif ForAll( elms, IsFFE ) then
> return FiniteFieldElements;
> elif ForAll( elms, IsPerm ) then
> return Permutations;
> elif ForAll( elms, IsMat ) then
> return Matrices;
> elif ForAll( elms, IsWord ) then
> return Words;
> elif ForAll( elms, IsAgWord ) then
> return AgWords;
> elif ForAll( elms, IsGaussInt ) then
> return GaussianIntegers;
> elif ForAll( elms, IsCyc ) then
> return Cyclotomics;
> else
> for elm in elms do
> if IsRec(elm) and IsBound(elm.domain)
> and ForAll( elms, l -> l in elm.domain )
> then
> return elm.domain;
> fi;
> od;
> Error("sorry, the elements lie in no common domain");
> fi;
> end;;
Of course we must define a function IsGaussInt, otherwise this could
not possibly work. This function is similar to the membership test we
already defined above.
gap> IsGaussInt := function ( x )
> return IsCycInt( x ) and (NofCyc( x ) = 1 or NofCyc( x ) = 4);
> end;;
Then we must define a function DefaultRing for the Gaussian integers
that does nothing but return GaussianIntegers.
gap> GaussianIntegersOps.DefaultRing := function ( elms )
> return GaussianIntegers;
> end;;
Now we can call Gcd with two Gaussian integers without having to pass
GaussianIntegers as first argument.
gap> Gcd( 2, 5 - E(4) );
1+E(4)
Of course GAP can not read your mind. In the following example it
assumes that you want to factor 10 over the ring of integers, not over
the ring of Gaussian integers (because Integers is the default ring
containing 10). So if you want to factor a rational integer over the
ring of Gaussian integers you must pass GaussianIntegers as first
argument.
gap> Factors( 10 );
[ 2, 5 ]
gap> Factors( GaussianIntegers, 10 );
[ -1-E(4), 1+E(4), 1+2*E(4), 2+E(4) ]
This concludes our example. In the file LIBNAME/"gaussian.g" you
will also find the definition of the field of Gaussian rationals. It is
so similar to the above definition that there is no point in discussing
it here. The next section shows you what further considerations are
necessary when implementing a type of parametrized domains (demonstrated
by implementing full symmetric permutation groups). For further details
see chapter Gaussians for a description of the Gaussian integers and
rationals and chapter Rings for a list of all functions applicable to
rings.
GAP 3.4.4