function main()
    print("Enter a size of square matrix: ");
    n = int(scan());
    a = [];
    println("Enter elements of matrix:");
    i = 0;
    while (i < n)
        # Enter i-th row of matrin
        append(a, []);
        j = 0;
        while (j < n)
            print("Enter a[", i, "][", j, "]: ");
            x = double(scan());
            # println("x = ", x);
            append(a[i], x);
            j = j + 1;
        endwhile
        i = i + 1;
    endwhile

    println("Matrix:");
    printMatrix(a);

    b = gauss(a);

    println("Row echelon form of matrix:");
    printMatrix(b);

    det = 1.;
    i = 0;
    while (i < n)
        det = det*b[i][i];
        i = i + 1;
    endwhile
    println("determinant = ", det);
endfunction

function printMatrix(a)
    m = size(a);
    i = 0;
    while (i < m)
        n = size(a[i]);
        j = 0;
        while (j < n)
            if (j > 0)
                print(" "); # separator
            endif
            print(a[i][j]);
            j = j + 1;
        endwhile
        println();
        i = i + 1;
    endwhile
endfunction

# Compute a row echelon form of a matrix
func gauss(a)

    eps = 1e-8;
    m = size(a);

    if (m == 0)
        return a;
    endif
    n = size(a[0]);
    if (n == 0)
        return a;
    endif

    i = 0;
    j = 0;
    while (i < m && j < n)
        # Find a maximal element in j-th column,
        #      starting from i-th row
        kmax = i; amax = fabs(a[kmax][j]);
        k = i + 1;
        while (k < m)
            if (fabs(a[k][j]) > amax)
                kmax = k;
                amax = fabs(a[k][j]);
            endif
            k = k + 1;
        endwhile
        if (amax <= eps)
            # Nearly zero column, erase very small numbers
            k = i;
            while (k < m)
                a[k][j] = 0.;
                k = k + 1;
            endwhile
            j = j + 1;  # Skip zero column
            continue;
        endif

        if (kmax != i)
            # Swap lines i and kmax
            s = j;
            while (s < n)
                tmp = a[i][s];
                a[i][s] = a[kmax][s];
                a[kmax][s] = -tmp;  # Minus to save determinant sign
                s = s + 1;
            endwhile
        endif

        # Annulate the j-th column, starting from row i+1
        k = i + 1;
        while (k < m)
            lambda = a[k][j]/a[i][j];
            # Subtract i-th row multiplied by lambda from k-th row
            a[k][j] = 0.;
            s = j + 1;
            while (s < n)
                a[k][s] = a[k][s] - lambda*a[i][s];
                s = s + 1;
            endwhile
            k = k + 1;
        endwhile

        # Go to the next minor of the matrix
        i = i + 1; 
        j = j + 1;
    endwhile
    return a;
endfunc

main();