chyanog 发表于 2019-11-23 18:21:58

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

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

#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;
}

mathe 发表于 2019-11-23 19:47:51

http://oeis.org/A000410

chyanog 发表于 2019-11-24 12:47:11

问题拓展一下，不限(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

chyanog 发表于 2019-12-3 22:12:12

本帖最后由 chyanog 于 2019-12-3 22:14 编辑

chyanog 发表于 2019-11-24 12:47
问题拓展一下，不限(0, 1)矩阵，比如(-1, 0, 1)，行列式换成积和式（Permanent），发现OIES已经收录了很多 ...

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

unsigned __int64 count;

int per(int **A, int n)
{
        if (n == 1)
                return A;
        if (n == 2)
                return (A * A + A * A);
        if (n == 3)
                return A * (A * A + A * A) + A * (A * A + A * A) + A * (A * A + A * A);
        if (n == 4)
                return A * (A * (A * A + A * A) + A * (A * A + A * A) + A * (A * A + A * A)) + A * (A * (A * A + A * A) + A * (A * A + A * A) + A * (A * A + A * A)) + A * (A * (A * A + A * A) + A * (A * A + A * A) + A * (A * A + A * A)) + A * (A * (A * A + A * A) + A * (A * A + A * A) + A * (A * A + A * A));
        int sum = 0, *m[--n];
        for (int i = 0; i < n; i++)
                m = A + 1;
        for (int i = 0; i <= n; i++) {
                sum += (A * per(m, n));
                if (i == n) break;
                m = A + 1;
        }
        return sum;
}

int permanent(int *A, int n)
{
        int *m;
        for (int i = 0; i < n; i++)
                m = A + (n * i);

        return per(m, n);
}

void f(int *a, int x, int n)
{
        for (a = -1; a <= 1; a++)
                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;
                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代码
n=4;
F=Symbol["i"<>IntegerString[#,10,2]]&;
perValue=Table,{i,n^2}]~Partition~n//Permanent;
iter=Table[{F,-1,1},{i,3,n^2}];

fun=Compile[{{i,_Integer}},
Module[{count=0},With[{i01=Floor-1},{i02=i-1-3Floor},Do,##2]];count],
RuntimeAttributes->{Listable},CompilationTarget->"C",RuntimeOptions->"Speed"
]&;

(len=Total@fun])//AbsoluteTiming

mathe 发表于 2019-12-4 18:17:43

两点改进意见:
i) 尽量避免递推调用，这个会有比较大的性能开销
ii) 交换矩阵任意两行或两列，不改变矩阵的奇异性。如果能够充分利用这个性质，可以大大降低计算复杂度。
页: [1]
查看完整版本: n×n的01奇异方阵的个数（代码优化）