// This is the Singular code that performs the computations described in Section 6
in the paper
//
// "Algebraic characterization of uniquely vertex colorable graphs"
//
// by Christopher J. Hillar and Troels Windfeldt.

// Options

option(redSB);
option(noredefine);
option(noprot);

timer = 1;

// Load libraries

LIB "general.lib";
LIB "qhmoduli.lib";

// Define functions

proc number_of_vertices(list G)
{
    return(G[1]);
}

proc number_of_edges(list G)
{
    return(nrows(G[2]));
}

proc edge_polynomial(list G, int n, int d)
{
    intmat edges = G[2];
    int i = edges[n,1];
    int j = edges[n,2];
    if (i < j)
    {
        return(x(i)^d-x(j)^d);
    }
    else
    {
        return(x(j)^d-x(i)^d);    
    }
}

proc graph_polynomial(list G)
{
    list l;
    int m = number_of_edges(G);
    for (int i = 1; i <= m; i++)
    {
        l[i] = edge_polynomial(G, i, 1);
    }
    return(product(l));
}

proc bin(int n, int k)
{
    list l = ringlist(basering);
    def R2 = ring(l);
    ring S = 0, x ,lp;
    int b = int(binomial(n, k, 0));
    setring(R2);
    return(b);
}

proc I(expr, int k)
{
    if (typeof(expr) == "int") { int n = expr; }
    else 
    {
        list G = expr;
        int n = number_of_vertices(G);
    }    
    ideal A;
    for (int i = 1; i <= n; i++)
    {
        A[i] = x(i)^k-1;
    }
    if (typeof(expr) == "int") { return(A); }
    int m = number_of_edges(G);
    for (int i = 1; i <= m; i++)
    {
        A[n+i] = quotient(edge_polynomial(G, i, k), edge_polynomial(G, i, 1))[1];
    }
    return(A);
}

proc J(int n, int k)
{
    if (k >= n) { return( ERROR("k must be less than n."); ); }
    k++;
    int m = bin(n,k);
    intmat B[m][k];
    for (int i = 1; i <= k; i++)
    {
        B[1,i] = i;    
    }
    for (int i = 2; i <= m; i++)
    {
        if (B[i-1,k] < n)
        {
            for (int j = 1; j < k; j++)
            {
                B[i,j] = B[i-1,j];
            }
            B[i,k] = B[i-1,k] + 1;
        }
        else
        {
            int l = k;
            while (B[i-1,l] == n-k+l) { l--; }
            for (int j = 1; j < l; j++) { B[i,j] = B[i-1,j]; }
            for (int j = l; j <= k; j++) { B[i,j] = B[i-1,l] + j - (l - 1); }
        }
    }
    ideal A;
    intmat edges[bin(k,2)][2];
    list G = n, edges;
    for (int l = 1; l <= nrows(B); l++)
    {
        int m = 1;
        for (int i = 1; i < k; i++)
        {
            for (int j = i + 1; j <= k; j++)
            {
                G[2][m,1] = B[l,i];
                G[2][m,2] = B[l,j];
                m++;
            }
        }
        A[l] = graph_polynomial(G);
    }
    return(A);
}

proc is_colorable(list G, int k, int m)
{
    if (m == 2)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        ideal A = I(n,k);
        ideal B = I(G,k);
        for (int i = 1; i <= m; i++)
        {
            A = ideal(groebner(A + ideal(B[n+i])));
        }
        A = std(A);
        return(vdim(A) != 0);    
    }
    if (m == 3)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        ideal A = I(n,k);
        ideal B = I(G,k);
        for (int i = 1; i <= m; i++)
        {
            A = ideal(groebner(A + ideal(B[n+i])));
        }
        A = std(A);
        return(reduce(1, A) != 0);
    }
    if (m == 4)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        poly f = 1;
        ideal A = I(n,k);
        for (int i = 1; i <= m; i++)
        {
            f = reduce(f * edge_polynomial(G, i, 1), A);
        }
        return(f != 0);
    }
    if (m == 5)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        poly f = 1;
        ideal A = J(n,k);
        for (int i = 1; i <= m; i++)
        {
            f = reduce(f * edge_polynomial(G, i, 1), A);
        }
        return(f != 0);
    }
    return( ERROR("Method m must be 2, 3, 4, or 5."); );
}

proc max_in_color_class(list G, int i)
{
    int n = number_of_vertices(G);
    list l;
    for (int j = 1; j <= n; j++)
    {
        if (G[3][j] == G[3][i]) { l = l + list(j); }
    }
    return(Max(l));
}

proc h(list U, int d)
{
    if (size(U) == 1) { return(x(U[1])^d); }
    else
    {
        return(quotient((h(delete(U, 2), d+1) - h(delete(U, 1),
d+1)),(x(U[1])-x(U[2])))[1]);
    }
}

proc g(list G, int i)
{
    int n = number_of_vertices(G);
    int k = Max(G[3]);
    if (i == n) { return(x(n)^k-1); }
    list l;
    for (int j = 1; j <= n; j++) { l = l + list(max_in_color_class(G, j)); }
    int m_0 = Min(l);
    int m_i = max_in_color_class(G, i);
    if (i == m_i or m_0 == m_i)
    {
        list v;
        list U;
        for (int j = 1; j <= n; j++) { v[j] = 0; }
        for (int j = 1; j <= n; j++)
        {
            int m_j = max_in_color_class(G, j);
            if (m_j > m_i) { v[m_j] = 1; }
        }
        v[i] = 1;
        for (int j = 1; j <= n; j++)
        {
            if (v[j] == 1) { U = U + list(j); }
        }
        return(h(U, k+1-size(U)));
    }
    else
    {
        return(x(i)-x(max_in_color_class(G, i)));
    }
}

proc g_bar(list G, int i)
{
    int n = number_of_vertices(G);
    int k = Max(G[3]);
    if (i == n) { return(1); }
    list l;
    for (int j = 1; j <= n; j++) { l = l + list(max_in_color_class(G, j)); }
    int m_0 = Min(l);
    int m_i = max_in_color_class(G, i);
    if (i == m_i or m_0 == m_i)
    {
        list v;
        list U;
        for (int j = 1; j <= n; j++) { v[j] = 0; }
        for (int j = 1; j <= n; j++)
        {
            int m_j = max_in_color_class(G, j);
            if (m_j > m_i) { v[m_j] = 1; }
        }
        v[i] = 1;
        for (int j = 1; j <= n; j++)
        {
            if (v[j] == 1) { U = U + list(j); }
        }
        intmat edges[bin(size(U),2)][2];
        list G = n, edges;
        int m = 1;
        for (int j1 = 1; j1 < size(U); j1++)
        {
            for (int j2 = j1 + 1; j2 <= size(U); j2++)
            {
                G[2][m,1] = U[j1];
                G[2][m,2] = U[j2];
                m++;
            }
        }
        return(graph_polynomial(G));
    }
    else
    {
        list U = i, max_in_color_class(G, i);
        return(h(U, 2));
    }
}

proc is_uniquely_colorable(list G, int m)
{
    int k = Max(G[3]);
    if (m == 3)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        ideal A = I(n,k);
        ideal B = I(G,k);
        for (int i = 1; i <= m; i++)
        {
            A = ideal(groebner(A + ideal(B[n+i])));
        }
        A = std(A);
        for (int i = 1; i <= n; i++)
        {
            if (reduce(g(G, i), A) != 0) { return(0); }
        }
        return(1);
    }
    if (m == 4)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        poly f = 1;
        poly l = 1;
        for (int i = 1; i <= n; i++) {l = l * g_bar(G, i); }
        ideal A = groebner(I(n,k) + ideal(l));
        for (int i = 1; i <= m; i++)
        {
            f = reduce(f * edge_polynomial(G, i, 1), A);
        }
        return(f == 0);        
    }
    if (m == 5)
    {
        int n = number_of_vertices(G);
        int m = number_of_edges(G);
        ideal A = I(n,k);
        ideal B = I(G,k);
        for (int i = 1; i <= m; i++)
        {
            A = ideal(groebner(A + ideal(B[n+i])));
        }
        A = std(A);
        return(string(vdim(A)) == factorial(k));    
    }
    return( ERROR("Method m must be 3, 4, or 5."); );
}

// Define the graph in Figure 1.

int n1 = 12;
int m1 = 23;
intmat edges1[m1][2] =
1,2,
1,4,
1,6,
1,12,
2,3,
2,5,
2,7,
3,8,
3,10,
4,9,
4,11,
5,6,
5,9,
5,12,
6,7,
6,10,
7,8,
7,11,
8,9,
8,12,
9,10,
10,11,
11,12;
intvec coloring1 = 1, 3, 2, 2, 1, 3, 2, 1, 3, 1, 3, 2;
list G1 = n1, edges1, coloring1;

// Define the graph in Figure 2.

int n2 = 24;
int m2 = 45;
intmat edges2[m2][2] =
1,2,
1,4,
1,6,
1,12,
2,3,
2,7,
3,10,
4,9,
4,11,
5,6,
5,9,
5,12,
6,7,
6,10,
7,8,
7,11,
8,9,
8,12,
10,11,
11,12,
13,14,
13,16,
13,18,
13,24,
14,15,
14,19,
15,22,
16,21,
16,23,
17,18,
17,21,
17,24,
18,19,
18,22,
19,20,
19,23,
20,21,
20,24,
22,23,
23,24,
2,14,
10,14,
10,22,
7,15,
3,18;
intvec coloring2 = 1, 3, 2, 2, 1, 3, 2, 1, 3, 1, 3, 2, 3, 2, 1, 1, 2, 1, 3, 2, 3, 3,
2, 1;
list G2 = n2, edges2, coloring2;

list G = G2;

// Main

ring R = 2, x(1..G[1]), dp;

system("--ticks-per-sec", "100");
int t = timer;

t = timer;
is_uniquely_colorable(G, 5);
t = timer-t; t;