Integral Determinant

Write a program to find the determinant of an integral square matrix. Note that the determinant of a square matrix can be defined recursively as follows: the determinant of a 11 matrix M = (a_1,1) is just the value |M| = a_1,1; further, the determinant of an matrix is

Here the notation M_1,i is the matrix by removing the first row and the i^th column of the original matrix M.

A straightforward method of calculating the determinant of an matrix by the recursive method will end up with n! multiplications, a very time-consuming algorithm. To give you a feeling about this, note that 15! = 1,307,674,368,000. To reduce the time complexity, there are two ways of modifying the original matrix for easier computation.

1.

Exchanging two columns (or rows) of a matrix will change the sign of the determinant; for example

2.

Multiplying one column by any constant, and add them to another column will not change the value of the determinant; for example:

Using the above methods, you shall be able to write a program for computing the determinants of matrices, even for a size like 3030, very efficiently. Below is an example to show how this can be done:

Note that the answer shall be an integer. That is, all the operations needed are just integer operations; by reducing to floating numbers would result in the round-off errors, which will be considered as the wrong answer. Do not worry about the problem of integral overflows problem. You can assume that the given data set will not cause the integer overflow problem. What is emphasized here is the required integer precision. Anyway use of floating numbers is not forbidden.

Input 

em Several sets of integral matrices. The inputs are just a list of integers. Within each set, the first integer (in a single line) represents the size of the matrix, n, which can be as large as 30, indication an matrix. After n, there will be n lines representing the n rows of the matrix; each line (row) contains exactly n integers. Thus, there is totally n² integers for the particular matrix. These matrices will occur repeatedly in the input as the pattern described above. An integer n = 0 (zero) signifies the end of input.

Output 

Sample Input 

2
5 2
3 4
3
2 3 5
1 6 7
4 8 9
0

Sample Output 

14
-27
*

計算行列式的值，題目已經先描述答案不會超過 int32 的範圍。
 不能使用降階法 O(n!)，因為 n 最大 30。
 一個解決方法，使用逆元做一次高斯消去，得到一個正整數。
 再使用一個浮點數的方法做一次高斯消去判斷原本行列式的正負，補回去即可。
 為了使用逆元，我們必須找到一個超過 int32 的大質數，但是乘起來又不能超過 int64(機率)。

#include <stdio.h>
#include <math.h>
#define LL long long
#define mod 2147483647LL
LL inverse(LL x, LL p) {
    if(!p)  return 1;
    if(p&1) return x*inverse((x*x)%mod, p>>1)%mod;
    return inverse((x*x)%mod, p>>1);
}
int doubleDet(int n, float A[][35]) {
    int inch = 0, i, j, k;
    float tmp;
    for(i = 0; i < n; i++) {
        int ch = -1;
        for(j = i; j < n; j++) {
            if(fabs(A[j][i]) >= 1e-6) {
                ch = j; break;
            }
        }
        if(ch == -1)
            return 0;
        if(i != ch) {
            for(j = i; j < n; j++)
                tmp = A[i][j], A[i][j] = A[ch][j], A[ch][j] = tmp;
            inch++;
        }
        LL inv, sub;
        for(j = i+1; j < n; j++) {
            if(fabs(A[j][i]) < 1e-6)  continue;
            tmp = A[j][i]/A[i][i];
            for(k = 0; k < n; k++) {
                A[j][k] = A[j][k] - A[i][k]*tmp;
            }
        }
    }
    tmp = 1;
    if(inch&1)  tmp *= -1;
    for(i = 0; i < n; i++)
        tmp *= A[i][i];
    return tmp > 0 ? 1 : -1;
}
int main() {
    int n, i, j, k;
    LL mix[35][35] = {}, tmp;
    float A[35][35] = {};
    while(scanf("%d", &n) == 1 && n) {
        for(i = 0; i < n; i++) {
            for(j = 0; j < n; j++) {
                scanf("%lld", &mix[i][j]);
                mix[i][j] = (mix[i][j]+mod)%mod;
                A[i][j] = mix[i][j];
            }
        }
        int inch = 0;
        for(i = 0; i < n; i++) {
            int ch = -1;
            for(j = i; j < n; j++) {
                if(mix[j][i]) {
                    ch = j; break;
                }
            }
            if(ch == -1 || !mix[ch][i]) {
                puts("0");  break;
            }
            if(i != ch) {
                for(j = i; j < n; j++)
                    tmp = mix[i][j], mix[i][j] = mix[ch][j], mix[ch][j] = tmp;
                inch++;
            }
            LL inv, sub;
            for(j = i+1; j < n; j++) {
                if(mix[j][i] == 0)  continue;
                inv = inverse(mix[i][i], mod-2);
                sub = (mix[j][i]*inv)%mod;
                for(k = 0; k < n; k++) {
                    mix[j][k] = mix[j][k] - sub*mix[i][k]%mod;
                    mix[j][k] = (mix[j][k]+mod)%mod;
                }
            }
        }
        if(i != n)  continue;
        LL ans = 1;
        for(i = 0; i < n; i++) {
            ans *= mix[i][i];
            ans %= mod;
        }
        if(inch&1)  ans *= -1;
        ans = (ans+mod)%mod;
        if(ans == 0 || doubleDet(n, A) >= 0)
            printf("%lld\n", ans);
        else
            printf("%lld\n", -(mod-ans));
    }
    puts("*");
    return 0;
}