数学研发论坛

 找回密码
 欢迎注册
楼主: 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<STATE;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, 2019-10-19 16:00 , Processed in 0.063262 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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