View Single Post
Old 2009-11-12, 07:29   #3
Dougy
 
Dougy's Avatar
 
Aug 2004
Melbourne, Australia

15210 Posts
Default

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 <stdio.h>
#include <gmp.h>
#include <time.h>

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.
Dougy is offline   Reply With Quote