找回密码
 欢迎注册
楼主: mathe

[讨论] 毒酒问题(加强版)

  [复制链接]
 楼主| 发表于 2009-6-10 09:42:53 | 显示全部楼层
这样判断过多了. 程序里面已经枚举了所有的a<=b,a<=c,c<=d的情况,而A,B,C,D枚举了所有的组合 也就是对于任何的(aA,bB),(cC,dD)组合已经穷举完了. 而如你那样修改,会在sa[a]|sa[b]) !=(sa[c]|sa[d])的情况下还要额外要求对应的la不等
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-6-10 10:18:00 | 显示全部楼层
sorry 没注意前面是 for(c=a;c
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-6-10 16:30:33 | 显示全部楼层
数学上的构造方法更需要天才的想法和可靠的理论
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-6-10 17:40:52 | 显示全部楼层
现在将这个模型的期望数目重新推导了一下,得出一个比较准确的公式,然后用计算机计算 的确在N=32时无法超过1000了(只有824),而N=33时期望数目可以达到1262. 所以这个方法33已经是最好的结果了.
  1. #define N 32
  2. #define MAXR 800
  3. double maxr;
  4. double f(double r[4],double x,int n)
  5. {
  6. double rt=0.0;
  7. int i;
  8. for(i=0;i<n;i++){
  9. rt*=x;
  10. rt+=r[n-1-i];
  11. }
  12. return rt;
  13. }
  14. double round(double x)
  15. {
  16. return (double)(int)(x+0.5);
  17. }
  18. #define ERROR 0.0001
  19. void solve3(double r[4],int m[])
  20. {
  21. double dr[3];
  22. double r0=100.0;
  23. double err;
  24. int c;
  25. int gc=0;
  26. dr[0]=r[1];dr[1]=2*r[2];dr[2]=3*r[3];
  27. start:
  28. c=0;
  29. do{
  30. double t=f(r,r0,4);
  31. err=t/f(dr,r0,3);
  32. r0-=err;
  33. c++;
  34. }while(fabs(err)>ERROR&&c<100);
  35. if(c>=100){
  36. r0=(double)rand();
  37. gc++;
  38. if(gc>=100){
  39. m[0]=m[1]=m[2]=-1;
  40. return;
  41. }
  42. goto start;
  43. }
  44. m[0]=(int)round(r0);
  45. dr[2]=r[3];
  46. dr[1]=r[2]-r0*dr[2];
  47. dr[0]=r[1]-r0*dr[1];
  48. double dd=dr[1]*dr[1]-4*dr[0]*dr[2];
  49. if(dd<0.0){
  50. m[1]=m[2]=-1;
  51. }else{
  52. dd=sqrt(dd);
  53. m[1]=(int)round((-dr[1]+dd)/(2*dr[2]));
  54. m[2]=(int)round((-dr[1]-dd)/(2*dr[2]));
  55. }
  56. }
  57. void get_m(int k,double q, int m[3])
  58. {
  59. int L=N-3*k;
  60. double r[4];
  61. double K=k*pow(3.0,(double)(k-2));
  62. double H=1-2*q*q+2*q*q*q;
  63. double Q=1-2*q*q+2*q*q*q*q;
  64. double HL=pow(H,L);
  65. double QL=pow(Q,L);
  66. double S=pow(6.0,k-2);
  67. double R=k*6*S+k*(k-1)/2.0*S*9.0;
  68. r[0]=-K;
  69. r[0]+=HL*K;
  70. r[1]=-3*HL*K;
  71. r[2]=3.0/2.0*HL*K;
  72. r[0]-=3.0/2.0*QL*K;
  73. r[1]+=11.0/2.0*QL*K;
  74. r[2]-=9.0/2.0*QL*K;
  75. r[3]=QL*K;
  76. r[1]+=0.5*QL*K*(K-1);
  77. r[2]-=3.0/2.0*QL*K*(K-1);
  78. r[3]+=QL*K*(K-1);
  79. r[3]+=QL*R*4.0;
  80. solve3(r,m);
  81. }
  82. void test(int k,double q, int m)
  83. {
  84. int L=N-3*k;
  85. double K=k*pow(3.0,(double)(k-2));
  86. double H=1-2*q*q+2*q*q*q;
  87. double Q=1-2*q*q+2*q*q*q*q;
  88. double HL=pow(H,L);
  89. double QL=pow(Q,L);
  90. double dm=m;
  91. double S=pow(6.0,k-2);
  92. double R=k*6*S+k*(k-1)/2.0*S*9.0;
  93. double r=K*dm-K*dm*(dm-1)*(dm-2)/2.0*HL-dm*(dm-1)*(dm-2)*(dm-3)/4.0*K*QL
  94. -K*(K-1)*dm*dm*(dm-1)*(dm-1)/4.0*QL-R*dm*dm*dm*dm*QL;
  95. if(r>=MAXR){
  96. printf("k=%d,q=%f,m=%d,r=%f\n",k,q,m,r);
  97. if(r>maxr)maxr=r;
  98. }
  99. }
  100. int _tmain(int argc, _TCHAR* argv[])
  101. {
  102. int k,s;
  103. double q;
  104. for(k=3;k<=8;k++){
  105. for(q=0.6667;q<0.707;q+=0.0001){
  106. int m[3];
  107. get_m(k,q,m);
  108. for(s=0;s<3;s++){
  109. if(m[s]>0&&m[s]<100000){
  110. test(k,q,m[s]);
  111. }
  112. }
  113. }
  114. }
  115. printf("maxr=%f\n",maxr);
  116. return 0;
  117. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-6-11 08:37:25 | 显示全部楼层
现在写了一个穷举的程序. 对于囚犯数目n不超过5,仅可以识别n+1坛酒(2坛毒酒) 对于n=6,可以识别8坛酒 00 03 0C 15 1A 26 29 30 对于n=7,可以识别10坛酒: 00 03 05 09 12 24 38 48 54 62 对于n=8,至少可以识别13坛酒: 00 03 0C 15 26 38 49 52 60 8A 90 A1 C4 对于n=9,至少可以识别16坛酒: 000 003 024 038 048 062 00D 016 051 081 098 0C4 10A 110 121 144 n=10,至少20坛酒 000 003 024 038 048 062 00D 016 051 098 0C4 10A 110 181 231 254 28A 2C1 30C 322 n=11,至少可以26坛酒 016 051 098 0C4 10A 110 181 231 254 28A 2C1 30C 322 000 003 00D 024 038 048 062 428 485 4C2 521 540 600 n=12,可以32坛 0C4 10A 110 081 144 221 288 301 431 494 505 602 628 740 000 00D 024 038 048 051 062 092 252 4C2 832 886 909 A04 A41 B02 C18 C80 n=13,可以39坛 010A 0110 0081 0121 0144 0205 0494 0612 0640 0720 0904 0A11 0A20 0B02 0C06 0C08 0C21 0E80 0000 0003 000D 0016 0024 0048 0051 02B0 03C0 08A8 101C 1022 10E0 1188 1242 1280 1401 1428 1444 1812 1901 n=14,至少47坛 0110 0081 0121 0144 0205 0494 0222 0460 0630 0704 0809 0842 0A11 0B02 0C10 0E80 10C0 1212 1301 1408 1421 1890 1940 1A00 0000 0003 000D 0016 0024 0038 0048 0051 0098 010A 1406 2046 2084 2214 22C0 2320 2402 2488 2908 2C20 3009 3010 3060
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-6-11 08:45:04 | 显示全部楼层
这下在Sloane网站中找到位置了: https://oeis.org/A054961 根据A054961,有 a(9) >= 17. Here is a candidate solution: {011010001 001001010 101000110 100000001 000010110 110010100 010000011 000110000 001001101 010000100 110001000 000101100 100100010 001000000 010100101 111100000 000011001}. - Dmitry Kamenetsky, Jul 17 2017 a(10) >= 22. Here is a candidate solution: {1010110000 1110000000 0010011000 0110100001 1000010101 0011000101 1001000000 0000000010 0001110100 0001001000 0010101010 1100010010 0000010001 0000101101 0100010100 0000011110 0100101000 0001100011 0000100100 0101000001 0101000110 1010000011}. - Dmitry Kamenetsky, Jul 17 2017 a(11) >= 31. Here is a candidate solution: {00101010010 01100101000 00001101100 00001010101 00100110100 00101001001 01000100101 00010101001 10001100001 11010000010 01001001010 00000001111 01111000000 00110010001 10000111000 11100000001 10010000101 00011110000 01000011001 00100100011 01010010100 01000110010 10110100000 10101000100 11001010000 00011000110 10000010011 10000100110 10100001010 00010011010 00000000100}. - Dmitry Kamenetsky, Jul 21 2017 a(12) >= 46. Here is a candidate solution: {011000011000 000000001000 000110100001 111000000010 110100000100 000010110100 010100010010 000011000011 010001001010 000111001000 010000101100 001001100010 100010000110 000001010110 001010101000 010010100010 001001001001 101000100100 000001100101 100011100000 000100011100 001000010011 010100001001 101100001000 100110010000 001110000010 100001001100 010101100000 110001000001 110010001000 011000100001 001100110000 001000001110 000010001101 110000110000 000100100110 000010011010 101001010000 000101010001 001100000101 100000010101 000001111000 101010000001 100101000010 010000000111 010010010001}. - Dmitry Kamenetsky, Jul 21 2017
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-6-14 17:26:04 | 显示全部楼层
请教mathe一个问题 对于n=6,可以识别8坛酒 00 03 0C 15 1A 26 29 30 列表中的00,13,...表示什么意思?? 即6个人怎么喝这8坛酒?/
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-6-14 21:38:21 | 显示全部楼层
上面每一行数都是一个16进制数,其中A,B,C,D,E,F表示16进制中的10~15. 而我们可以将每位16进制数转化为4位二进制数,比如 0 -> 0000 1 -> 0001 2 -> 0010 3 -> 0011 4 -> 0100 ... F -> 1111 于是替换以后,每一行可以表示成一个n位二进制数.其中每一位分别对应那个人是否要喝对应的那坛酒(所在行表示酒) 如上面n=6的情况,那么第一行00,表示第一坛酒6个人都不喝.第二行0C转化为6位二进制数是001100,表示第3和第4个人喝第二坛酒.第三行15即二进制010101,所以第三坛酒是第2,4,6个人喝,等等.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-6-15 08:34:16 | 显示全部楼层
百度贴巴中的tannis_jin关于这个题目找到了一种挺好的构造方法(虽然这个方法对于大规模数据的结果还是不够好),这里我来整理一下他的思路: 假设h(m)表示m坛酒中如果有1坛或2坛毒酒的情况下,需要囚犯的最少数目. 那么对于一个m*m方阵坛酒, 首先,我们可以将每一行的酒混合,得到m坛混合酒,其中有1坛或2坛有毒,所以需要h(m)个囚犯可以确定在哪一航. 同样在使用h(m)个囚犯可以确定在哪一列. 在确定行和列以外,在任意一条对角线方向,除了对角线顶点,以及紧挨着顶点的酒坛,余下的酒按桶对角线垂直方向的直线分类可以分成2m-5类,我们最多使用h(2m-5)个囚犯可以确定在那条对角线上. 于是采用这个方案有 $h(m*m)<=2h(m)+h(2m-5)$ 而对于1000坛酒,我们可以取m=32,于是$h(1024)<=2h(32)+h(59)$ 我们看65#的结果,其中每个结果中都包含一个全0的结果,去掉那个全0的方案就正好对应h(.)的方案 也就是说包含1到2坛毒酒的38坛只需要13个人识别,所以h(32)<=13.而看上面序列的趋势,h(59)应该不会超过16,也就是用这个方法需要2*13+16=42个人就可以了. 而如果我们继续结合7#的方案,那么一个边长为19的六边形里面包含1027桶酒.其中三组对边方向,每个方向有37条直线,所以我们只需要3h(37)=3*13=39个人酒可以确定了. 而将方法推广到m非常大的情况, tannis_jin的方法相当于使用递推式 $h(m*m)<=2h(m)+h(2m-5)$ 于是这个方法得到的h(m)的复杂度也为$h(m)=O(log^2(m))$ 而结合7#的方案当中,其递推式为 $h(3*m*(m-1)+1)=3h(m)$, 这个方案的复杂度为$h(m)=O(log^{log_2 3}(m))$,比上面的方法要好一些. 但是无论比起24#的集合方法(其中1000坛酒可以得到41人的结果),47#4n个人解决$2^n$坛酒的$F_{2^n}$上的构造方法(1000坛酒40人的结果),还是39#半随机方法(1000坛酒33人的结果59#)(它们的复杂度都是$O(log(n))$)都还是要差一些.可能这是构造方法本身的缺陷了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-6-15 11:33:59 | 显示全部楼层
mathe有相当的数学专业精神,佩服至极.....
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-11-23 18:07 , Processed in 0.024461 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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