/* Copyright (C) 2002 J. M. Spivey */

/** Solve the LightsOut puzzle */
public class LightsSolver {
    /* This solves a specific system of 25 simultaneous equations over Z3.  The
    equations arise from the Lights Out puzzle; the input is the vector y,
    and the output is a vector x that solves the equation M x = y, where M
    is a fixed matrix determined by the structure of the puzzle.

    Actually, the matrix M is singular, so there may be no solutions or
    several.  If one or more solutions exist, the solution with the
    smallest total is chosen. */

    public static final int N = 5, K = N*N;
    
    /** mod3 -- carefully reduce modulo 3 */
    private static int mod3(int i) {
        i %= 3;
        if (i < 0) i += 3;
        return i;
    }

    /** makeMatrix -- set up the matrix of coefficients */
    private static int[][] makeMatrix() {
        int M[][] = new int[K][K];
        
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++) {
                int k = N*i+j;
                M[k][k] = 1;
                if (i > 0) M[k][k-N] = 1;
                if (i < N-1) M[k][k+N] = 1;
                if (j > 0) M[k][k-1] = 1;
                if (j < N-1) M[k][k+1] = 1;
            }
        
        return M;
    }
    
    /** swap -- swap two elements of a vector */
    private static void swap(int v[], int i, int j) {
        int tmp = v[i]; v[i] = v[j]; v[j] = tmp;
    }

    /** swapCols -- swap two columns of a matrix */
    private static void swapCols(int M[][], int i, int j) {
        for (int k = 0; k < K; k++) {
            int tmp = M[k][i]; M[k][i] = M[k][j]; M[k][j] = tmp;
        }
    }
    
    /** swapRows -- swap two rows of a matrix */
    private static void swapRows(int M[][], int i, int j) {
        for (int k = 0; k < K; k++) {
            int tmp = M[i][k]; M[i][k] = M[j][k]; M[j][k] = tmp;
        }
    }
            
    /** gauss -- Gaussian elimination */
    private static int gauss(int M[][], int perm[], int y[]) {
        /* This reduces M to an upper-triangular matrix U such that
         LUP = QM, where Q and P are permutation matrices and L is
         lower-triangular.  Actually, U has the form
    
            ( A B )
            ( 0 0 )

         where A is unitary upper triangular r x r, and gauss returns 
         the rank r of M.  It also replaces y by a vector y' such that 
	 Ly' = Qy, and returns the permutation P as perm. 

         Subsequently, we can solve the equations by back substitution, as
         UPx = y' ==> QMx = LUPx = Ly' = Qy ==> Mx = y. */
        
        /* Initialize perm to the identity permutation */
        for (int p = 0; p < K; p++) perm[p] = p;

        for (int p = 0; p < K; p++) {
            /* Find a row r >= p and column c >= p with M[r][c] != 0 */
            int c = p, r = p;
            while (c < K && M[r][c] == 0) {
                r++;
                if (r >= K) { c++; r = p; }
            }
            if (c >= K) return p; // Singular

            /* Swap rows r and p, columns c and p */
            swapRows(M, r, p); swap(y, r, p);
            swapCols(M, c, p); swap(perm, c, p);

            /* Normalize row p by dividing (same as multiplying 
	       mod 3) by M[p][p] */
            int u = M[p][p]; M[p][p] = 1;
            for (int q = p+1; q < K; q++) M[p][q] = mod3(M[p][q] * u);
            y[p] = mod3(y[p] * u);

            /* Eliminate column p from rows [p+1..K) */
            for (r = p+1; r < K; r++) {
                int v = M[r][p]; 
                if (v == 0) continue;
                M[r][p] = 0;
                for (int q = p+1; q < K; q++) 
                    M[r][q] = mod3(M[r][q] - v * M[p][q]);
                y[r] = mod3(y[r] - v * y[p]);
            }
        }
        
        return K;
    }

    /** Derive unknown values in x by back-substitution */
    private static void backSubst(int M[][], int x[], 
				  int perm[], int y[], int rank) {
        for (int p = rank-1; p >= 0; p--) {
            int u = y[p];
            for (int q = p+1; q < K; q++)
                u = mod3(u - M[p][q] * x[perm[q]]);
            x[perm[p]] = u;
        }
    }

    /** weight -- sum of components of a vector */
    private static int weight(int v[]) {
        int z = 0;
        for (int i = 0; i < K; i++) z += v[i];
        return z;
    }

    /** solve -- Find smallest solution */
    public static int[] solve(int y[]) {
        int M[][] = makeMatrix();
        int perm[] = new int[K];
        int rank = gauss(M, perm, y);
        
        for (int i = rank; i < K; i++)
            if (y[i] != 0) return null; // No solution
        
        // Compute number of solutions n
        int n = 1;
        for (int i = rank; i < K; i++) n *= 3;
        
        // Try each solution, keep the one with least weight
        int best[] = null;
        int bweight = 1000000;
        
        for (int s = 0; s < n; s++) {
            int x[] = new int[K];
            
            // Compute the free variables
            int k = s;
            for (int i = K-1; i >= rank; i--) {
                x[perm[i]] = k % 3;
                k /= 3;
            }
            
            // Derive other x values by back substitution
            backSubst(M, x, perm, y, rank);

            // Keep this solution if it is the best so far
            int z = weight(x);
            if (best == null || z < bweight) {
                best = x; bweight = z;
            }
        }
        
        return best;
    }
}