//////////////////////////////////////////////////////////////////////
//		Magma routines for diptych varieties
//
// Gavin Brown and Miles Reid
// August 2012
//
// References are to:
// Diptych Varieties, I.
//
//////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////
// Compute the matrices (r a // b s) and (r g // h s) from d, e, k
// arising in Theorem 3.3.
//     * d,e,k as in the theorem.
//     * type := 1 for the first factorisation, equation (3.5)
//       type := 2 for the second, equation (3.6).
//////////////////////////////////////////////////////////////////////

function DiptychMatrixPair(d,e,k:type:=1)
    D1 := Matrix(2,2,[Integers()|d,-1,1,0]);
    E1 := Matrix(2,2,[Integers()|e,-1,1,0]);
    D2 := Matrix(2,2,[Integers()|0,-1,1,d]);
    E2 := Matrix(2,2,[Integers()|0,-1,1,e]);
    twist := Matrix(2,2,[Integers()|0,1,-1,0]);
    M := IdentityMatrix(Integers(),2);
    N := IdentityMatrix(Integers(),2);
    for i in [1..k] do
        if IsOdd(i) then
            M := M*D1;
            N := E2*N;
        else
            M := M*E1;
            N := D2*N;
        end if;
    end for;
    M := M*twist;
    N := twist*N;
    if type eq 1 then
        return M,N;
    else
        return N,M;
    end if;
end function;


//////////////////////////////////////////////////////////////////////
// Compute various sequences of tags given the matrix pair M,N of
// the Classification Theorem.
// The tags are computed in Section 2.
//////////////////////////////////////////////////////////////////////

// The continued fraction coefficients (with minus signs) of rational q.

function DiptychContinuedFraction(q)
    Q := [Integers()|];
    while not IsIntegral(q) do
	n := Ceiling(q);
    	Append(~Q,n);
	q := 1/(n-q);
    end while;
    Append(~Q,Integers()!q);
    return Q;
end function;

// The continued fraction expansion as a rational number
//	a0 - 1/(a1 - 1/(a2 - ... )) for A = [a0,a1,...].

function DiptychExpansion(A)
    if #A eq 1 then
    	return A[1];
    else
        return A[1] - 1/$$(A[2..#A]);
    end if;
end function;

// The tags on the long rectangle from a diptych pair of matrices
// M and N that come from Theorem 3.3.

function DiptychSequencesFromMatrixPair(M,N)
    r := M[1,1];
    a := M[1,2];
    b := M[2,1];
    s := M[2,2];
    g := N[1,2];
    h := N[2,1];
    // Check that this is a matrix pair
    conditions := (N[1,1] eq r) and (N[2,2] eq s) and
    	((a+h) mod r eq 0) and ((a+h) mod s eq 0) and
    	((b+g) mod r eq 0) and ((b+g) mod s eq 0);
    error if not conditions, "This is not a matrix pair";
    d := (a+h) div s;
    e := (b+g) div r;
    printf "d = %o\te = %o\n",d,e;
    printf "a0 = %o\tak = %o\n",Ceiling(-b/r),Ceiling(a/r);
    print "-----";
    printf "[a0..ak-1] = -b/r = %o\n[ak..a1] = a/r = %o\n",
        DiptychContinuedFraction(-b/r), DiptychContinuedFraction(a/r);
    printf "\t[ak-1..a1] = r/(r-al) = %o\n",
        DiptychContinuedFraction(r/(r - (a mod r)));
    printf "[b0..bl-1] = -a/s = %o\n[bl..b1] = b/s = %o\n",
        DiptychContinuedFraction(-a/s), DiptychContinuedFraction(b/s);
    printf "\t[bl-1..b1] = s/(s-be) = %o\n",
        DiptychContinuedFraction(s/(s - (b mod s)));
    print "-----";
    // printf "[a0'..ak-1] = -b/r = %o\n[ak'..a1] = a/r = %o\n",
    printf "[a0'..ak-1] = g/r = %o\n[ak'..a1] = -h/r = %o\n",
        DiptychContinuedFraction(g/r), DiptychContinuedFraction(-h/r);
    printf "[b0'..bl-1] = h/s = %o\n[bl'..b1] = -g/s = %o\n",
        DiptychContinuedFraction(h/s), DiptychContinuedFraction(-g/s);
    print "-----";
    return "";
end function;


//////////////////////////////////////////////////////////////////////
/*		Testing

// Extended example from the paper:
M := Matrix(2,2,[7,12,4,7]);
N := Matrix(2,2,[7,24,2,7]);
DiptychSequencesFromMatrixPair(M,N);

// An example of the main case:
> d := 2;
> e := 4;
> k := 4;
> M,N := DiptychMatrixPair(d,e,k);
> M,N;                            
[12 41]
[ 7 24]

[12 41]
[ 7 24]
> DiptychSequencesFromMatrixPair(M,N);
d = 2   e = 4
a0 = 0  ak = 4
-----
[a0..ak-1] = -b/r = [ 0, 2, 4, 2 ]
[ak..a1] = a/r = [ 4, 2, 4, 2 ]
        [ak-1..a1] = r/(r-al) = [ 2, 4, 2 ]
[b0..bl-1] = -a/s = [ -1, 2, 2, 4, 2, 2 ]
[bl..b1] = b/s = [ 1, 2, 2, 4, 2, 2 ]
        [bl-1..b1] = s/(s-be) = [ 2, 2, 4, 2, 2 ]
-----
[a0'..ak-1] = g/r = [ 4, 2, 4, 2 ]
[ak'..a1] = -h/r = [ 0, 2, 4, 2 ]
[b0'..bl-1] = h/s = [ 1, 2, 2, 4, 2, 2 ]
[bl'..b1] = -g/s = [ -1, 2, 2, 4, 2, 2 ]
-----


> d := 2;
> e := 4;
> k := 5;
> M,N := DiptychMatrixPair(d,e,k);
> M,N;                            
[41 70]
[24 41]

[ 41 140]
[ 12  41]
> DiptychSequencesFromMatrixPair(M,N);
d = 2   e = 4
a0 = 0  ak = 2
-----
[a0..ak-1] = -b/r = [ 0, 2, 4, 2, 4 ]
[ak..a1] = a/r = [ 2, 4, 2, 4, 2 ]
        [ak-1..a1] = r/(r-al) = [ 4, 2, 4, 2 ]
[b0..bl-1] = -a/s = [ -1, 2, 2, 4, 2, 3 ]
[bl..b1] = b/s = [ 1, 3, 2, 4, 2, 2 ]
        [bl-1..b1] = s/(s-be) = [ 3, 2, 4, 2, 2 ]
-----
[a0'..ak-1] = g/r = [ 4, 2, 4, 2, 4 ]
[ak'..a1] = -h/r = [ 0, 4, 2, 4, 2 ]
[b0'..bl-1] = h/s = [ 1, 2, 2, 4, 2, 3 ]
[bl'..b1] = -g/s = [ -3, 3, 2, 4, 2, 2 ]
-----

// Same example but with matrices switched - so long rectangles swap
// and turn upside down:
>   d := 2;
>   e := 4;
>   k := 5;
>   M,N := DiptychMatrixPair(d,e,k: type:=2);
> M,N;                                       
[ 41 140]
[ 12  41]

[41 70]
[24 41]
>   DiptychSequencesFromMatrixPair(M,N);
d = 4   e = 2
a0 = 0  ak = 4
-----
[a0..ak-1] = -b/r = [ 0, 4, 2, 4, 2 ]
[ak..a1] = a/r = [ 4, 2, 4, 2, 4 ]
        [ak-1..a1] = r/(r-al) = [ 2, 4, 2, 4 ]
[b0..bl-1] = -a/s = [ -3, 3, 2, 4, 2, 2 ]
[bl..b1] = b/s = [ 1, 2, 2, 4, 2, 3 ]
        [bl-1..b1] = s/(s-be) = [ 2, 2, 4, 2, 3 ]
-----
[a0'..ak-1] = g/r = [ 2, 4, 2, 4, 2 ]
[ak'..a1] = -h/r = [ 0, 2, 4, 2, 4 ]
[b0'..bl-1] = h/s = [ 1, 3, 2, 4, 2, 2 ]
[bl'..b1] = -g/s = [ -1, 2, 2, 4, 2, 3 ]
-----


// A pair of matrices with e=1, so the alternating game is not visible:
> d := 5;
> e := 1;
> k := 6;
> M,N := DiptychMatrixPair(d,e,k);
> M,N;
[40 29]
[11  8]

[40 29]
[11  8]
> DiptychSequencesFromMatrixPair(M,N);
d = 5   e = 1
a0 = 0  ak = 1
-----
[a0..ak-1] = -b/r = [ 0, 4, 3, 4 ]
[ak..a1] = a/r = [ 1, 4, 3, 4 ]
        [ak-1..a1] = r/(r-al) = [ 4, 3, 4 ]
[b0..bl-1] = -a/s = [ -3, 2, 3, 2 ]
[bl..b1] = b/s = [ 2, 2, 3, 2 ]
        [bl-1..b1] = s/(s-be) = [ 2, 3, 2 ]
-----
[a0'..ak-1] = g/r = [ 1, 4, 3, 4 ]
[ak'..a1] = -h/r = [ 0, 4, 3, 4 ]
[b0'..bl-1] = h/s = [ 2, 2, 3, 2 ]
[bl'..b1] = -g/s = [ -3, 2, 3, 2 ]
-----

*/
//////////////////////////////////////////////////////////////////////

