1.28 About Defining New Domains

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.

Previous Up Top Next
Index

GAP 3.4.4
April 1997