数学研发论坛

 找回密码
 欢迎注册
查看: 186|回复: 6

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

[复制链接]
发表于 2019-11-23 18:21:58 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

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

  1. #include <stdio.h>
  2. #include <time.h>

  3. int main()
  4. {
  5.         clock_t t0 = clock();
  6.         int count = 0;
  7.         for (int i = 0; i < 1<<25; i++)
  8.         {
  9.                 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;
  10.                 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;
  11.                 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)));
  12.                 if (det==0)
  13.                         count++;
  14.         }
  15.         printf("%d\n",count);
  16.         printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);
  17.         return 0;
  18. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-11-23 19:47:51 来自手机 | 显示全部楼层
http://oeis.org/A000410
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 3 天前 | 显示全部楼层
本帖最后由 chyanog 于 2019-12-3 22:14 编辑
chyanog 发表于 2019-11-24 12:47
问题拓展一下,不限(0, 1)矩阵,比如(-1, 0, 1),行列式换成积和式(Permanent),发现OIES已经收录了很多 ...


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

  3. unsigned __int64 count;

  4. int per(int **A, int n)
  5. {
  6.         if (n == 1)
  7.                 return A[0][0];
  8.         if (n == 2)
  9.                 return (A[0][0] * A[1][1] + A[0][1] * A[1][0]);
  10.         if (n == 3)
  11.                 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]);
  12.         if (n == 4)
  13.                 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]));
  14.         int sum = 0, *m[--n];
  15.         for (int i = 0; i < n; i++)
  16.                 m[i] = A[i + 1] + 1;
  17.         for (int i = 0; i <= n; i++) {
  18.                 sum += (A[i][0] * per(m, n));
  19.                 if (i == n) break;
  20.                 m[i] = A[i] + 1;
  21.         }
  22.         return sum;
  23. }

  24. int permanent(int *A, int n)
  25. {
  26.         int *m[n];
  27.         for (int i = 0; i < n; i++)
  28.                 m[i] = A + (n * i);

  29.         return per(m, n);
  30. }

  31. void f(int *a, int x, int n)
  32. {
  33.         for (a[n] = -1; a[n] <= 1; a[n]++)
  34.                 if (n == 0) {
  35.                         if (permanent(a, x) == 0)
  36.                                 count++;
  37.                 }
  38.                 else
  39.                         f(a, x, n - 1);
  40. }

  41. int main()
  42. {
  43.         clock_t t0 = clock();
  44.         for (int i = 1; i <= 4; i++) {
  45.                 count = 0;
  46.                 int a[i*i];
  47.                 f(a, i, i*i - 1);
  48.                 printf("%llu\n", count);
  49.         }

  50.         printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);
  51.         return 0;
  52. }

  53. // 1, 33, 7555, 13482049, 186481694371
复制代码

Mathematica代码
  1. n=4;
  2. F=Symbol["i"<>IntegerString[#,10,2]]&;
  3. perValue=Table[F[i],{i,n^2}]~Partition~n//Permanent;
  4. iter=Table[{F[i],-1,1},{i,3,n^2}];

  5. fun=Compile[{{i,_Integer}},
  6. Module[{count=0},With[{i01=Floor[i/3]-1},{i02=i-1-3Floor[i/3]},Do[If[#==0,count++],##2]];count],
  7. RuntimeAttributes->{Listable},CompilationTarget->"C",RuntimeOptions->"Speed"
  8. ]&[perValue,Sequence@@iter];

  9. (len=Total@fun[Range[0,8]])//AbsoluteTiming
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 前天 18:17 | 显示全部楼层
两点改进意见:
i) 尽量避免递推调用,这个会有比较大的性能开销
ii) 交换矩阵任意两行或两列,不改变矩阵的奇异性。如果能够充分利用这个性质,可以大大降低计算复杂度。

点评

是的,要避免递归  发表于 昨天 12:58
谢谢。递推应该是递归吧。  发表于 昨天 11:23
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-12-6 18:03 , Processed in 0.130385 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表