n×n的01奇异方阵的个数（代码优化）

公子相

http://oeis.org/A046747
统计行列式为0的方阵即可，计算5×5及以上的速度就开始慢了。下面的代码是计算5阶的，GCC需要1秒左右，VC 2秒左右。应该怎么优化？除了对称性，生成所有01的组合怎么做更高效，计算行列式使用一个显式的长表达式速度应该不比循环慢吧

1. #include <stdio.h>#include <time.h>int main() {        clock_t t0 = clock();        int count = 0;        for (int i = 0; i < 1<<25; i++)        {                int i01, i02, i03, i04, i05, i06, i07, i08, i09, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, i21, i22, i23, i24, i25;                i01 = (i >> 0) & 1; i02 = (i >> 1) & 1; i03 = (i >> 2) & 1; i04 = (i >> 3) & 1; i05 = (i >> 4) & 1; i06 = (i >> 5) & 1; i07 = (i >> 6) & 1; i08 = (i >> 7) & 1; i09 = (i >> 8) & 1; i10 = (i >> 9) & 1; i11 = (i >> 10) & 1; i12 = (i >> 11) & 1; i13 = (i >> 12) & 1; i14 = (i >> 13) & 1; i15 = (i >> 14) & 1; i16 = (i >> 15) & 1; i17 = (i >> 16) & 1; i18 = (i >> 17) & 1; i19 = (i >> 18) & 1; i20 = (i >> 19) & 1; i21 = (i >> 20) & 1; i22 = (i >> 21) & 1; i23 = (i >> 22) & 1; i24 = (i >> 23) & 1; i25 = (i >> 24) & 1;                int det = i05*(-i09*(i13*(i16*i22-i17*i21)-i12*(i16*i23-i18*i21)+i11*(i17*i23-i18*i22))+i08*(i14*(i16*i22-i17*i21)-i12*(i16*i24-i19*i21)+i11*(i17*i24-i19*i22))-i07*(i14*(i16*i23-i18*i21)-i13*(i16*i24-i19*i21)+i11*(i18*i24-i19*i23))+i06*(i14*(i17*i23-i18*i22)-i13*(i17*i24-i19*i22)+i12*(i18*i24-i19*i23)))-i04*(-i10*(i13*(i16*i22-i17*i21)-i12*(i16*i23-i18*i21)+i11*(i17*i23-i18*i22))+i08*(i15*(i16*i22-i17*i21)-i12*(i16*i25-i20*i21)+i11*(i17*i25-i20*i22))-i07*(i15*(i16*i23-i18*i21)-i13*(i16*i25-i20*i21)+i11*(i18*i25-i20*i23))+i06*(i15*(i17*i23-i18*i22)-i13*(i17*i25-i20*i22)+i12*(i18*i25-i20*i23)))+i03*(-i10*(i14*(i16*i22-i17*i21)-i12*(i16*i24-i19*i21)+i11*(i17*i24-i19*i22))+i09*(i15*(i16*i22-i17*i21)-i12*(i16*i25-i20*i21)+i11*(i17*i25-i20*i22))-i07*(i15*(i16*i24-i19*i21)-i14*(i16*i25-i20*i21)+i11*(i19*i25-i20*i24))+i06*(i15*(i17*i24-i19*i22)-i14*(i17*i25-i20*i22)+i12*(i19*i25-i20*i24)))-i02*(-i10*(i14*(i16*i23-i18*i21)-i13*(i16*i24-i19*i21)+i11*(i18*i24-i19*i23))+i09*(i15*(i16*i23-i18*i21)-i13*(i16*i25-i20*i21)+i11*(i18*i25-i20*i23))-i08*(i15*(i16*i24-i19*i21)-i14*(i16*i25-i20*i21)+i11*(i19*i25-i20*i24))+i06*(i15*(i18*i24-i19*i23)-i14*(i18*i25-i20*i23)+i13*(i19*i25-i20*i24)))+i01*(-i10*(i14*(i17*i23-i18*i22)-i13*(i17*i24-i19*i22)+i12*(i18*i24-i19*i23))+i09*(i15*(i17*i23-i18*i22)-i13*(i17*i25-i20*i22)+i12*(i18*i25-i20*i23))-i08*(i15*(i17*i24-i19*i22)-i14*(i17*i25-i20*i22)+i12*(i19*i25-i20*i24))+i07*(i15*(i18*i24-i19*i23)-i14*(i18*i25-i20*i23)+i13*(i19*i25-i20*i24)));                if (det==0)                        count++;        }        printf("%d\n",count);        printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);        return 0;}

复制代码

epezzaqegazeu

http://oeis.org/A000410

ucobiit

问题拓展一下，不限(0, 1)矩阵，比如(-1, 0, 1)，行列式换成积和式（Permanent），发现OIES已经收录了很多相关的数列，但是还没有收录 Number of n X n (-1, 0,1)-matrices with zero permanent. 好像可以提交一个新数列了。
前5项为：
1, 33, 7555, 13482049, 186481694371

litong

已提交该数列，可以由以下C代码计算。计算积和式目前还没有比较高效的算法。经测试展开后的表达式确实比用循环快（相当于循环展开），但是大于四阶时式子太长，就不再展开了。下面的C代码还不如Mathematica版速度快。

1. #include <stdio.h>#include <time.h>unsigned __int64 count;int per(int **A, int n){        if (n == 1)                return A[0][0];        if (n == 2)                return (A[0][0] * A[1][1] + A[0][1] * A[1][0]);        if (n == 3)                return A[0][0] * (A[1][1] * A[2][2] + A[1][2] * A[2][1]) + A[0][1] * (A[1][0] * A[2][2] + A[1][2] * A[2][0]) + A[0][2] * (A[1][0] * A[2][1] + A[1][1] * A[2][0]);        if (n == 4)                return A[0][0] * (A[1][1] * (A[2][2] * A[3][3] + A[2][3] * A[3][2]) + A[1][2] * (A[2][1] * A[3][3] + A[2][3] * A[3][1]) + A[1][3] * (A[2][1] * A[3][2] + A[2][2] * A[3][1])) + A[0][1] * (A[1][0] * (A[2][2] * A[3][3] + A[2][3] * A[3][2]) + A[1][2] * (A[2][0] * A[3][3] + A[2][3] * A[3][0]) + A[1][3] * (A[2][0] * A[3][2] + A[2][2] * A[3][0])) + A[0][2] * (A[1][0] * (A[2][1] * A[3][3] + A[2][3] * A[3][1]) + A[1][1] * (A[2][0] * A[3][3] + A[2][3] * A[3][0]) + A[1][3] * (A[2][0] * A[3][1] + A[2][1] * A[3][0])) + A[0][3] * (A[1][0] * (A[2][1] * A[3][2] + A[2][2] * A[3][1]) + A[1][1] * (A[2][0] * A[3][2] + A[2][2] * A[3][0]) + A[1][2] * (A[2][0] * A[3][1] + A[2][1] * A[3][0]));        int sum = 0, *m[--n];        for (int i = 0; i < n; i++)                m[i] = A[i + 1] + 1;        for (int i = 0; i <= n; i++) {                sum += (A[i][0] * per(m, n));                if (i == n) break;                m[i] = A[i] + 1;        }        return sum;}int permanent(int *A, int n){        int *m[n];        for (int i = 0; i < n; i++)                m[i] = A + (n * i);        return per(m, n);}void f(int *a, int x, int n){        for (a[n] = -1; a[n] <= 1; a[n]++)                if (n == 0) {                        if (permanent(a, x) == 0)                                count++;                }                else                        f(a, x, n - 1);}int main(){        clock_t t0 = clock();        for (int i = 1; i <= 4; i++) {                count = 0;                int a[i*i];                f(a, i, i*i - 1);                printf("%llu\n", count);        }        printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);        return 0;}// 1, 33, 7555, 13482049, 186481694371

复制代码

Mathematica代码

1. n=4;F=Symbol["i"<>IntegerString[#,10,2]]&;perValue=Table[F[i],{i,n^2}]~Partition~n//Permanent;iter=Table[{F[i],-1,1},{i,3,n^2}];fun=Compile[{{i,_Integer}},Module[{count=0},With[{i01=Floor[i/3]-1},{i02=i-1-3Floor[i/3]},Do[If[#==0,count++],##2]];count],RuntimeAttributes->{Listable},CompilationTarget->"C",RuntimeOptions->"Speed"]&[perValue,Sequence@@iter];(len=Total@fun[Range[0,8]])//AbsoluteTiming

复制代码

xiaobai

两点改进意见:
i) 尽量避免递推调用，这个会有比较大的性能开销
ii) 交换矩阵任意两行或两列，不改变矩阵的奇异性。如果能够充分利用这个性质，可以大大降低计算复杂度。