The algorithms explained in the paper have been implemented in Magma and can be downloaded here.
Note that currently all elements and subgroups of the modular groups have to be given as elements and subgroups of SL(2, Z). The following example shows the main features of our implementation.
###
How to use this file with Magma:

Attach the downloaded file to your running Magma process.
```
Attach("sl2z.m");
```

To check if some subgroup H of PSL(2,Z) has finite index in PSL(2,Z), one uses the IndexInPSL2Z intrinsic.
```
> H:= sub< SL(2, Integers()) | [1,4,0,1], [0,1,-1,0] >;
> IndexInPSL2Z(H);
Infinity
> G:= sub< SL(2, Integers()) | [1,2,0,1], [1,0,2,1] >;
> IndexInPSL2Z(G);
6
```

So the group H has infinite index in the modular group while G has index 6. To decide if G is a congruence subgroup, one can use MyIsCongruence.
```
> MyIsCongruence(G);
true
```

So G is indeed a congruence subgroup.
To test whether an element of PSL(2,Z) lies in H one can use the Membership command.
```
> M:= SL(2, Integers()) ! [1,2,0,1];
> ok:= Membership(H, M); ok;
false
```

If one turns the CosetRep flag on, the second return value will be some data representing the right coset. This data can be turned into a matrix with CosetToMatrix. The coset data as well as the matrix are unique and can therefore be used for comparing cosets.
```
> M:= SL(2, Integers()) ! [1,2,0,1];
> ok, coset:= Membership(H, M : CosetRep);
> CosetToMatrix(coset);
[ 1 2]
[-4 -7]
> ok, coset2:= Membership(H, H.1 * M : CosetRep); coset eq coset2;
true
```

Suppose H is given by n generators and h in H.
If the Word flag is turned on, the Membership intrinsic will return a straight line progam (SLP) representing h as the third return value.
Evaluating this SLP in H itself yields (not surprisingly) h. Evaluating the SLP in a free group of rank n expresses h as a word in the n generators of H.
This solves the constructive membership problem.
Here is an example.
```
> M:= SL(2, Integers()) ! [ -40, -159, -1, -4 ];
> ok, _, w:= Membership(H, M : Word);
> assert ok and Evaluate(w, H) in {M, -M};
> Evaluate( w, FreeGroup(Ngens(H)) );
$.1^2 * $.2^2 * $.1^8 * $.2^-1 * $.1
> H.1^2 * H.2^2 * H.1^8 * H.2^-1 * H.1 in {M, -M};
true
```