> 
# These procedures implement basic arithmetic for Eisenstein integers 
# represented as pairs of rational integers. The complex number a+b*w, 
# where w:=(1+sqrt(-3))/2, is represented as [a,b]. 
# 
# Long division of Eisenstein integers can be done in different ways. 
# Here, I used an algorithm that doesn't bother to take the nearest 
# Eisenstein integer, but simply rounds to the nearest integer in the pair 
# representation. This procedure may well cause violations of algebraic 
# rules one would expect, such as that long division of a into b should 
# give the same result as long division of c*a into c*b, whenever c is a 
# unit, and should give the same quotient for any c. There is a lot of 
# symmetry, and it complicates these details. 
> eisentimes:=proc(u,v)local a,b,c,d;\
a:=u[1];b:=u[2];c:=v[1];\
d:=v[2];[a*c-b*d,b*c+(a+b)*d]end;

eisentimes :=

    proc(u,v)
    local a,b,c,d;
        a := u[1]; b := u[2]; c := v[1]; d := v[2]; [a*c-b*d,b*c+(a+b)*d]
    end

> eisenplus:=proc(u,v)local a,b,c,d;\
a:=u[1];b:=u[2];c:=v[1];d:=v[2];[a+c,b+d]end;

eisenplus := proc(u,v)
             local a,b,c,d;
                 a := u[1]; b := u[2]; c := v[1]; d := v[2]; [a+c,b+d]
             end

> eisenminus:=proc(u,v)local a,b,c,d;\
a:=u[1];b:=u[2];c:=v[1];d:=v[2];[a-c,b-d]end;

eisenminus := proc(u,v)
              local a,b,c,d;
                  a := u[1]; b := u[2]; c := v[1]; d := v[2]; [a-c,b-d]
              end

> eisencon:=proc(u)[u[1]+u[2],-u[2]]end;\
#complex conjugate;

eisencon := proc(u) [u[1]+u[2],-u[2]] end

> 
> eisentimes([0,1],[0,1]);

                                     [-1, 1]
> eisendiv:=proc(u,v)local a,b,c,d,e,f,g,h;\
a:=u[1];b:=u[2];c:=v[1];d:=v[2];\
e:=eisentimes(u,eisencon(v));\
f:=eisentimes(v,eisencon(v))[1];\
g:=floor(e[1]/f+1/2);\
h:=floor(e[2]/f+1/2);\
[[g,h],eisenplus(u,eisentimes([-g,-h],v))]end;

eisendiv := proc(u,v)
            local a,b,c,d,e,f,g,h;
                a := u[1];
                b := u[2];
                c := v[1];
                d := v[2];
                e := eisentimes(u,eisencon(v));
                f := eisentimes(v,eisencon(v))[1];
                g := floor(e[1]/f+1/2);
                h := floor(e[2]/f+1/2);
                [[g,h],eisenplus(u,eisentimes([-g,-h],v))]
            end

> 
> eisendiv([7,5],[2,1]);

                                [[4, 0], [-1, 1]]
> eisentimes([4,0],[2,1]);

                                     [8, 4]
> dot:=proc(u,v)\
eisenplus(eisentimes(u[1],v[1]),eisentimes(u[2],v[2]))end;

dot := proc(u,v) eisenplus(eisentimes(u[1],v[1]),eisentimes(u[2],v[2])) end

> 
# This is for later use in keeping track of gcd calculations using a 
# matrix. One way to think of gcd calculations is that you're building a 
# 2 by 2 matrix of the form [[a,b],[c,d]] so that all along it has 
# determinant 1, and at the end, the bottom row is the original numbers 
# u,v  with their gcd taken out, and the top row is  numbers n m so that 
# m*u-n*v=+-1. At each stage you multiply the old matrix by the 
# simple matrix [[q,1],[1,0]] where q is the quotient in the simple gcd 
# algorithm.  Of course, with our representation of Eisenstein integers 
# as pairs, we can't just call up the library routines for matrix 
# multiplication in the E-context. The routine immediately below 
# implements this idea for the rational integers. 
> gcdmatshowntell:=proc(u,v) local m,uu,vv,q,t;\
m:=array([[1,0],[0,1]]);uu:=u;vv:=v;\
while vv<>0 do \
printf(`Now u,v, and the matrix are\n`);\
print([uu,vv,[[m[1,1],m[1,2]],[m[2,1],m[2,2]]]]);\
q:=floor(uu/vv);t:=uu mod vv;uu:=vv;vv:=t;\
m:=multiply(array([[q,1],[1,0]]),m)od;\
[uu,vv,[[m[1,1],m[1,2]],[m[2,1],m[2,2]]]]end;

gcdmatshowntell := proc(u,v)
                   local m,uu,vv,q,t;
                       m := array([[1,0],[0,1]]);
                       uu := u;
                       vv := v;
                       while vv <> 0 do
                           printf(`Now u,v, and the matrix are
                                  `                           );
                           print([uu,vv,[[m[1,1],m[1,2]],[m[2,1],m[2,2]]]]);
                           q := floor(uu/vv);
                           t := uu mod vv;
                           uu := vv;
                           vv := t;
                           m := multiply(array([[q,1],[1,0]]),m)
                       od;
                       [uu,vv,[[m[1,1],m[1,2]],[m[2,1],m[2,2]]]]
                   end

> with(linalg):gcdmatshowntell(24,17);
Now u,v, and the matrix are


                           [24, 17, [[1, 0], [0, 1]]]
Now u,v, and the matrix are


                            [17, 7, [[1, 1], [1, 0]]]
Now u,v, and the matrix are


                            [7, 3, [[3, 2], [1, 1]]]
Now u,v, and the matrix are


                            [3, 1, [[7, 5], [3, 2]]]

                           [1, 0, [[24, 17], [7, 5]]]
> 
> dot([[1,1],[2,3]],[[1,1],[2,3]]);

                                    [-5, 24]
> 
> w:=(1+sqrt(-3))/2;

                                                1/2
                              w := 1/2 + 1/2 I 3
> (1+w)^2+(2+3*w)^2;

                                  1/2 2                 1/2 2
                    (3/2 + 1/2 I 3   )  + (7/2 + 3/2 I 3   )
> expand(");

                                            1/2
                                  7 + 12 I 3
> -5+24*w;

                                            1/2
                                  7 + 12 I 3
> eisentimes([1,1],[1,1]);

                                     [0, 3]
> expand((1+w)^2);

                                             1/2
                                3/2 + 3/2 I 3
> 3*w;

                                             1/2
                                3/2 + 3/2 I 3
> eisentimes([2,3],[2,3]);

                                    [-5, 21]
> expand((2+3*w)^2);

                                              1/2
                               11/2 + 21/2 I 3
> -5+21*w;

                                              1/2
                               11/2 + 21/2 I 3
> matx:=proc(u,v)local\
u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4;\
u1:=u[1][1];u2:=u[1][2];u3:=u[2][1];\
u4:=u[2][2];v1:=v[1][1];v2:=v[1][2];v3:=v[2][1];v4:=v[2][2];\
w1:=dot([u1,u2],[v1,v3]);\
w2:=dot([u1,u3],[v2,v4]);w3:=dot([u3,u4],[v1,v3]);\
w4:=dot([u3,u4],[v2,v4]);[[w1,w2],[w3,w4]]end;\
#matx implements mult'n of 2 by 2 Eisenstein matrices

matx := proc(u,v)
        local u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4;
            u1 := u[1][1];
            u2 := u[1][2];
            u3 := u[2][1];
            u4 := u[2][2];
            v1 := v[1][1];
            v2 := v[1][2];
            v3 := v[2][1];
            v4 := v[2][2];
            w1 := dot([u1,u2],[v1,v3]);
            w2 := dot([u1,u3],[v2,v4]);
            w3 := dot([u3,u4],[v1,v3]);
            w4 := dot([u3,u4],[v2,v4]);
            [[w1,w2],[w3,w4]]
        end

> 
> r:=[[1,2],[3,4]];

                              r := [[1, 2], [3, 4]]
> r[2,2];
Error, invalid subscript selector

> 
# Note: In MapleVR4 they let you access lists-of-lists entries in this 
# fashion. That is, the same two lines won't give you an error message. 
> r[2][2];

                                        4
# Inelegant, but it works. 
> one:=[1,0];zero:=[0,0];

                                  one := [1, 0]

                                 zero := [0, 0]
> id:=[[one,zero],[zero,one]];

                   id := [[[1, 0], [0, 0]], [[0, 0], [1, 0]]]
> f1:=[[[2,3],one],[one,zero]];

                   f1 := [[[2, 3], [1, 0]], [[1, 0], [0, 0]]]
> matx(f1,id);

                      [[[2, 3], [1, 0]], [[1, 0], [0, 0]]]
> matx(f1,f1);

                     [[[-4, 21], [2, 3]], [[2, 3], [1, 0]]]
> matx(f1,f1);

                     [[[-4, 21], [2, 3]], [[2, 3], [1, 0]]]
> matx(f1,");

                   [[[-69, 96], [-4, 21]], [[-4, 21], [2, 3]]]
> matx(",");

          [[[-4880, -3759], [-1811, 276]], [[-1811, 276], [-430, 294]]]
> eisengcd:=proc(u,v)local uu,vv,qu,rem,a,b,m;\
uu:=u;vv:=v;m:=id;\
while vv<>zero do \
a:=eisendiv(uu,vv);qu:=a[1];rem:=a[2];\
m:=matx(m,[[qu,one],[one,zero]]);\
uu:=vv;vv:=rem\
od;[uu,m]end;

eisengcd := proc(u,v)
            local uu,vv,qu,rem,a,b,m;
                uu := u;
                vv := v;
                m := id;
                while vv <> zero do
                    a := eisendiv(uu,vv);
                    qu := a[1];
                    rem := a[2];
                    m := matx(m,[[qu,one],[one,zero]]);
                    uu := vv;
                    vv := rem
                od;
                [uu,m]
            end


eisengcd := proc(u,v)
            local uu,vv,qu,rem,a,b,m;
                uu := u;
                vv := v;
                m := id;
                while vv <> zero do
                    a := eisendiv(uu,vv);
                    qu := a[1];
                    rem := a[2];
                    m := matx(m,[[qu,one],[one,zero]]);
                    uu := vv;
                    vv := rem
                od;
                [uu,m]
            end

> eg:=eisengcd([7,5],[3,2]);

            eg := [[0, -1], [[[-12, 7], [7, 0]], [[-5, 3], [3, 0]]]]
> eisennorm:=proc(u) eisentimes(u,eisencon(u))[1]end;

eisennorm := proc(u) eisentimes(u,eisencon(u))[1] end

> eisennorm([7,5]);

                                       109
> eisennorm([-12,7]);

                                       109
> eisendet:=proc(m)local m11,m12,m21,m22;m11:=m[1][1];m12:=m[1][2];m21:=m[2][1];m22:=m[2][2];eisenminus(eisentimes(m11,m22),eisentimes(m21,m12))end;

eisendet := proc(m)
            local m11,m12,m21,m22;
                m11 := m[1][1];
                m12 := m[1][2];
                m21 := m[2][1];
                m22 := m[2][2];
                eisenminus(eisentimes(m11,m22),eisentimes(m21,m12))
            end


edet :=

    proc(m) eisenplus(eisentimes(m[1,1],m[2,2]),-eisentimes(m[1,2],m[2,1])) end

> eg2:=eg[2];

                 eg2 := [[[-12, 7], [7, 0]], [[-5, 3], [3, 0]]]
> eisendet(eg2);

                                     [-1, 0]
> eisendiv([7,5],[-12,7]);

                                [[0, -1], [0, 0]]
> eisengcd([107,58],[231,99]);

        [[0, -1], [[[-165, 107], [81, -60]], [[-330, 231], [161, -128]]]]
> eisendiv([107,58],[-165,107]);

                                [[0, -1], [0, 0]]
> 
>