Thread: Counting Latin rectangles View Single Post
 2009-11-12, 07:29 #3 Dougy     Aug 2004 Melbourne, Australia 15210 Posts Thanks (:. They say the thesis will be available online, so I can post a link. Well not much of a race to respond... I'm not sure if it's because it's quite off-topic for this forum. The code, however, should be able to be improved without knowledge of Latin rectangles. Perhaps I'll try to explain how it works a bit better. I'm mostly interested in R(5,n) since for R(4,n) I have lots of data already and R(6,n) seems too hard. Overall, it's a sum over "partitions" of n=A[0]+A[1]+...+A[15] (where the A[i] can take the values of 0...n). For each partition of n the variable "internal" is added to "count." Code: #include #include #include int n,x,y; time_t start_time; mpz_t count,internal,temp; int A[16],B[16]; int DoyleF5_1(int A[16]) { return A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]; } int DoyleF5_12(int A[16]) { return A[1]+A[2]+A[3]+A[4]; } int DoyleF5_123(int A[16]) { return A[1]+A[2]; } int DoyleF5_1234(int A[16]) { return A[1]; } int DoyleF5_124(int A[16]) { return A[1]+A[3]; } int DoyleF5_13(int A[16]) { return A[1]+A[2]+A[5]+A[6]; } int DoyleF5_134(int A[16]) { return A[1]+A[5]; } int DoyleF5_14(int A[16]) { return A[1]+A[3]+A[5]+A[7]; } int DoyleF5_2(int A[16]) { return A[1]+A[2]+A[3]+A[4]+A[9]+A[10]+A[11]+A[12]; } int DoyleF5_23(int A[16]) { return A[1]+A[2]+A[9]+A[10]; } int DoyleF5_234(int A[16]) { return A[1]+A[9]; } int DoyleF5_24(int A[16]) { return A[1]+A[3]+A[9]+A[11]; } int DoyleF5_3(int A[16]) { return A[1]+A[2]+A[5]+A[6]+A[9]+A[10]+A[13]+A[14]; } int DoyleF5_34(int A[16]) { return A[1]+A[5]+A[9]+A[13]; } int DoyleF5_4(int A[16]) { return A[1]+A[3]+A[5]+A[7]+A[9]+A[11]+A[13]+A[15]; } After the front matter, the above is required for Doyle's formula. In each case, it's just addition of some collection of the A[i]'s. Code: int DoyleG5(int A[16]) { return DoyleF5_1(A)*DoyleF5_2(A)*DoyleF5_3(A)*DoyleF5_4(A)-DoyleF5_1(A)*DoyleF5_2(A)*DoyleF5_34(A)-DoyleF5_1(A)*DoyleF5_23(A)*DoyleF5_4(A)+2*DoyleF5_1(A)*DoyleF5_234(A)-DoyleF5_1(A)*DoyleF5_24(A)*DoyleF5_3(A)-DoyleF5_12(A)*DoyleF5_3(A)*DoyleF5_4(A)+DoyleF5_12(A)*DoyleF5_34(A)+2*DoyleF5_123(A)*DoyleF5_4(A)-6*DoyleF5_1234(A)+2*DoyleF5_124(A)*DoyleF5_3(A)-DoyleF5_13(A)*DoyleF5_2(A)*DoyleF5_4(A)+DoyleF5_13(A)*DoyleF5_24(A)+2*DoyleF5_134(A)*DoyleF5_2(A)-DoyleF5_14(A)*DoyleF5_2(A)*DoyleF5_3(A)+DoyleF5_14(A)*DoyleF5_23(A); } This is the guts. But it's just multiplication, addition and subtraction of the earlier functions. Code: int MinusOneOrNot(int A[16]) { if((A[2]+A[3]+A[5]+A[8]+A[9]+A[12]+A[14]+A[15])%2==0) { return 1; } else { return -1; } } This will later decide the sign of "internal." Code:  start_time=time(NULL); mpz_init(count); mpz_init(temp); mpz_init(internal); for(n=20;n<=100;n++) { mpz_set_si(count,0); Initialisation... I've set it to run for n=20 to 100. Code:  for(A[1]=0;A[1]<=n;A[1]++) { for(A[2]=0;A[2]<=n-A[1];A[2]++) { for(A[3]=0;A[3]<=n-A[1]-A[2];A[3]++) { for(A[4]=0;A[4]<=n-A[1]-A[2]-A[3];A[4]++) { for(A[5]=0;A[5]<=n-A[1]-A[2]-A[3]-A[4];A[5]++) { for(A[6]=0;A[6]<=n-A[1]-A[2]-A[3]-A[4]-A[5];A[6]++) { for(A[7]=0;A[7]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6];A[7]++) { for(A[8]=0;A[8]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7];A[8]++) { for(A[9]=0;A[9]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8];A[9]++) { for(A[10]=0;A[10]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9];A[10]++) { for(A[11]=0;A[11]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10];A[11]++) { for(A[12]=0;A[12]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10]-A[11];A[12]++) { for(A[13]=0;A[13]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10]-A[11]-A[12];A[13]++) { for(A[14]=0;A[14]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10]-A[11]-A[12]-A[13];A[14]++) { for(A[15]=0;A[15]<=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10]-A[11]-A[12]-A[13]-A[14];A[15]++) { A[0]=n-A[1]-A[2]-A[3]-A[4]-A[5]-A[6]-A[7]-A[8]-A[9]-A[10]-A[11]-A[12]-A[13]-A[14]-A[15]; Chooses the partition of n=A[0]+A[1]+...+A[15]. Is there a better way? Code:  mpz_fac_ui(internal,n); for(x=0;x<=15;x++) { mpz_fac_ui(temp,A[x]); mpz_div(internal,internal,temp); } mpz_mul_si(internal,internal,MinusOneOrNot(A)); Initialising "internal" to be the multinomial coefficient ${n \choose A[0],A[1],\dots,A[15]}=n!/(A[0]! \cdot A[1]! \cdots A[15]!)$. Then MinusOneOrNot(A) determines the appropriate sign. Code:  for(x=0;x<=15;x++) { for(y=0;y<=15;y++) { B[y]=A[y]; } B[x]--; mpz_set_si(temp,0); mpz_set_si(temp,DoyleG5(B)); mpz_pow_ui(temp,temp,A[x]); mpz_mul(internal,internal,temp); } This is the product required for Doyle's formula. for each 0<=x<=15 the B[i]=A[i] except B[x] is reduced by one. "internal" is then created via a "temp" variable which is $\prod_{0 \leq x \leq 15} DoyleG5(B)^{A[x]}$ as a result of this loop. Code:  } mpz_add(count,count,internal); } } } } } } } } } } } } } } } if(n>=5) { mpz_div_ui(count,count,(n-1)*(n-2)*(n-3)*(n-4)); } gmp_printf("R(5,%d)=%Zd %ld sec.\n",n,count,time(NULL)-start_time); } } Adds "internal" to the variable "count" for each partition. At the end it divides by (n-1)(n-2)(n-3)(n-4) when it can and prints out the result. Hope this helps. But if there's anything I could explain more clearly, I'm happy to.