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

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

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[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代码

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

复制代码

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

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

账号		自动登录	找回密码
密码			欢迎注册

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

点评