Optimal a-set pair code

From Polymath Wiki
Jump to navigationJump to search

// Optimally selecting symmetric representatives of subsets of [2]^4

  1. include <stdio.h>
  2. include <stdlib.h>
  3. include <string.h>
  4. include <time.h>


  1. define NUM_PERMS 384

int perms[NUM_PERMS][16]; // The 48 symmetries of the 4D cube long pow2[27]; // The 27 powers of 2

void buildperms(void) { 

  int i,j,k,l,p,x,y,z,w,a,b,c,d;
  int n=0;
  long power=1;

  for (i=0; i<27; i++) { pow2[i] = power; power *= 2; }

  for (x=0;x<2;x++)

for (y=0;y<2;y++) for (z=0;z<2;z++) 

 		      for (w=0;w<2;w++)

for (p=0;p<24;p++) { for (i=0;i<2;i++) for (j=0;j<2;j++) for (k=0;k<2;k++) 

                         for (l=0;l<2;l++)

{ switch (p) { case 0: a=i; b=j; c=k; d=l; break; case 1: a=i; b=j; c=l; d=k; break; case 2: a=i; b=k; c=j; d=l; break; case 3: a=i; b=k; c=l; d=j; break; case 4: a=i; b=l; c=k; d=j; break; case 5: a=i; b=l; c=j; d=k; break;

case 6: a=j; b=i; c=k; d=l; break; case 7: a=j; b=i; c=l; d=k; break; case 8: a=j; b=k; c=i; d=l; break; case 9: a=j; b=k; c=l; d=i; break; case 10: a=j; b=l; c=k; d=i; break; case 11: a=j; b=l; c=i; d=k; break;

case 12: a=k; b=j; c=i; d=l; break; case 13: a=k; b=j; c=l; d=i; break; case 14: a=k; b=i; c=l; d=j; break; case 15: a=k; b=i; c=j; d=l; break; case 16: a=k; b=l; c=i; d=j; break; case 17: a=k; b=l; c=j; d=i; break;

case 18: a=l; b=j; c=k; d=i; break; case 19: a=l; b=j; c=i; d=k; break; case 20: a=l; b=k; c=j; d=i; break; case 21: a=l; b=k; c=i; d=j; break; case 22: a=l; b=i; c=k; d=j; break; case 23: a=l; b=i; c=j; d=k; break; } if (x) a = 1-a; if (y) b = 1-b; if (z) c = 1-c; if (w) d = 1-d; perms[n][i+2*j+4*k+8*l] = a+2*b+4*c+8*d; } n++; } }


long perm(long set, int i) // permutes set by permutation i { long in,out; int j;

in = set; out = 0; for (j=0; j < 16; j++) { if (in & 1) out |= pow2[perms[i][j]]; in >>= 1; }

return out; }


long astatnum[256];

void init_astats(void) // precomputed 

{

astatnum[0]=82173l; 

   astatnum[1]=82173l;
   astatnum[2]=82173l;
   astatnum[3]=51060l;
   astatnum[4]=82173l;
   astatnum[5]=51060l;
   astatnum[6]=51052l;
   astatnum[7]=16836l;
   astatnum[8]=82173l;
   astatnum[9]=51052l;
   astatnum[10]=51060l;
   astatnum[11]=16836l;
   astatnum[12]=51060l;
   astatnum[13]=16836l;
   astatnum[14]=16836l;
   astatnum[15]=5750l;
   astatnum[16]=82173l;
   astatnum[17]=51060l;
   astatnum[18]=51052l;
   astatnum[19]=16836l;
   astatnum[20]=51052l;
   astatnum[21]=16836l;
   astatnum[22]=20472l;
   astatnum[23]=3071l;
   astatnum[24]=69868l;
   astatnum[25]=25777l;
   astatnum[26]=25777l;
   astatnum[27]=4021l;
   astatnum[28]=25777l;
   astatnum[29]=4021l;
   astatnum[30]=5000l;
   astatnum[31]=625l;
   astatnum[32]=82173l;
   astatnum[33]=51052l;
   astatnum[34]=51060l;
   astatnum[35]=16836l;
   astatnum[36]=69868l;
   astatnum[37]=25777l;
   astatnum[38]=25777l;
   astatnum[39]=4021l;
   astatnum[40]=51052l;
   astatnum[41]=20472l;
   astatnum[42]=16836l;
   astatnum[43]=3071ll;
   astatnum[44]=25777l;
   astatnum[45]=5000l;
   astatnum[46]=4021l;
   astatnum[47]=625l;
   astatnum[48]=51060l;
   astatnum[49]=16836l;
   astatnum[50]=16836l;
   astatnum[51]=5750l;
   astatnum[52]=25777l;
   astatnum[53]=4021l;
   astatnum[54]=5000l;
   astatnum[55]=625l;
   astatnum[56]=25777l;
   astatnum[57]=5000l;
   astatnum[58]=4021l;
   astatnum[59]=625l;
   astatnum[60]=9604l;
   astatnum[61]=784l;
   astatnum[62]=784l;
   astatnum[63]=98l;
   astatnum[64]=82173l;
   astatnum[65]=51052l;
   astatnum[66]=69868l;
   astatnum[67]=25777l;
   astatnum[68]=51060l;
   astatnum[69]=16836l;
   astatnum[70]=25777l;
   astatnum[71]=4021l;
   astatnum[72]=51052l;
   astatnum[73]=20472l;
   astatnum[74]=25777l;
   astatnum[75]=5000l;
   astatnum[76]=16836l;
   astatnum[77]=3071l;
   astatnum[78]=4021l;
   astatnum[79]=625l;
   astatnum[80]=51060l;
   astatnum[81]=16836l;
   astatnum[82]=25777l;
   astatnum[83]=4021l;
   astatnum[84]=16836l;
   astatnum[85]=5750l;
   astatnum[86]=5000l;
   astatnum[87]=625l;
   astatnum[88]=25777l;
   astatnum[89]=5000l;
   astatnum[90]=9604l;
   astatnum[91]=784l;
   astatnum[92]=4021l;
   astatnum[93]=625l;
   astatnum[94]=784l;
   astatnum[95]=98l;
   astatnum[96]=51052l;
   astatnum[97]=20472l;
   astatnum[98]=25777l;
   astatnum[99]=5000l;
   astatnum[100]=25777l;
   astatnum[101]=5000l;
   astatnum[102]=9604l;
   astatnum[103]=784l;
   astatnum[104]=20472l;
   astatnum[105]=4825l;
   astatnum[106]=5000l;
   astatnum[107]=512l;
   astatnum[108]=5000l;
   astatnum[109]=512l;
   astatnum[110]=784l;
   astatnum[111]=64l;
   astatnum[112]=16836l;
   astatnum[113]=3071l;
   astatnum[114]=4021l;
   astatnum[115]=625l;
   astatnum[116]=4021l;
   astatnum[117]=625l;
   astatnum[118]=784l;
   astatnum[119]=98l;
   astatnum[120]=5000l;
   astatnum[121]=512l;
   astatnum[122]=784l;
   astatnum[123]=64l;
   astatnum[124]=784l;
   astatnum[125]=64l;
   astatnum[126]=64l;
   astatnum[127]=8l;
   astatnum[128]=82173l;
   astatnum[129]=69868l;
   astatnum[130]=51052l;
   astatnum[131]=25777l;
   astatnum[132]=51052l;
   astatnum[133]=25777l;
   astatnum[134]=20472l;
   astatnum[135]=5000l;
   astatnum[136]=51060l;
   astatnum[137]=25777l;
   astatnum[138]=16836l;
   astatnum[139]=4021l;
   astatnum[140]=16836l;
   astatnum[141]=4021l;
   astatnum[142]=3071l;
   astatnum[143]=625l;
   astatnum[144]=51052l;
   astatnum[145]=25777l;
   astatnum[146]=20472l;
   astatnum[147]=5000l;
   astatnum[148]=20472l;
   astatnum[149]=5000l;
   astatnum[150]=4825l;
   astatnum[151]=512l;
   astatnum[152]=25777l;
   astatnum[153]=9604l;
   astatnum[154]=5000l;
   astatnum[155]=784l;
   astatnum[156]=5000l;
   astatnum[157]=784l;
   astatnum[158]=512l;
   astatnum[159]=64l;
   astatnum[160]=51060l;
   astatnum[161]=25777l;
   astatnum[162]=16836l;
   astatnum[163]=4021l;
   astatnum[164]=25777l;
   astatnum[165]=9604l;
   astatnum[166]=5000l;
   astatnum[167]=784l;
   astatnum[168]=16836l;
   astatnum[169]=5000l;
   astatnum[170]=5750l;
   astatnum[171]=625l;
   astatnum[172]=4021l;
   astatnum[173]=784l;
   astatnum[174]=625l;
   astatnum[175]=98l;
   astatnum[176]=16836l;
   astatnum[177]=4021l;
   astatnum[178]=3071l;
   astatnum[179]=625l;
   astatnum[180]=5000l;
   astatnum[181]=784l;
   astatnum[182]=512l;
   astatnum[183]=64l;
   astatnum[184]=4021l;
   astatnum[185]=784l;
   astatnum[186]=625l;
   astatnum[187]=98l;
   astatnum[188]=784l;
   astatnum[189]=64l;
   astatnum[190]=64l;
   astatnum[191]=8l;
   astatnum[192]=51060l;
   astatnum[193]=25777l;
   astatnum[194]=25777l;
   astatnum[195]=9604l;
   astatnum[196]=16836l;
   astatnum[197]=4021l;
   astatnum[198]=5000l;
   astatnum[199]=784l;
   astatnum[200]=16836l;
   astatnum[201]=5000l;
   astatnum[202]=4021l;
   astatnum[203]=784l;
   astatnum[204]=5750l;
   astatnum[205]=625l;
   astatnum[206]=625l;
   astatnum[207]=98l;
   astatnum[208]=16836l;
   astatnum[209]=4021l;
   astatnum[210]=5000l;
   astatnum[211]=784l;
   astatnum[212]=3071l;
   astatnum[213]=625l;
   astatnum[214]=512l;
   astatnum[215]=64l;
   astatnum[216]=4021l;
   astatnum[217]=784l;
   astatnum[218]=784l;
   astatnum[219]=64l;
   astatnum[220]=625l;
   astatnum[221]=98l;
   astatnum[222]=64l;
   astatnum[223]=8l;
   astatnum[224]=16836l;
   astatnum[225]=5000l;
   astatnum[226]=4021l;
   astatnum[227]=784ll;
   astatnum[228]=4021l;
   astatnum[229]=784l;
   astatnum[230]=784l;
   astatnum[231]=64l;
   astatnum[232]=3071l;
   astatnum[233]=512l;
   astatnum[234]=625l;
   astatnum[235]=64l;
   astatnum[236]=625l;
   astatnum[237]=64l;
   astatnum[238]=98l;
   astatnum[239]=8l;
   astatnum[240]=5750l;
   astatnum[241]=625l;
   astatnum[242]=625l;
   astatnum[243]=98l;
   astatnum[244]=625l;
   astatnum[245]=98l;
   astatnum[246]=64l;
   astatnum[247]=8l;
   astatnum[248]=625l;
   astatnum[249]=64l;
   astatnum[250]=98l;
   astatnum[251]=8l;
   astatnum[252]=98l;
   astatnum[253]=8l;
   astatnum[254]=8l;
   astatnum[255]=1l;
}


// const int a[16] = {0, 2, 6, 8, 18, 20, 24, 26}; // coeffs of a-points

int main() { 

  int i,j,k,l,m;

  buildperms();
  init_astats();

  unsigned long set, reflected, bestset, deflect;
  long double count,best;
  long a, b;
  unsigned int mult;
  int numa;

  long double cost = 0, cost2=0;
  long double num_apairs = 0;



  for (set = 0; set < pow2[16]; set++)

{ numa = 0; for (i=0; i<16; i++) if (set & pow2[i]) numa++; // a = number of a-points if (numa==1 || numa==2 || numa==3) continue;

mult=0; a = astatnum[set/256]; b = astatnum[set%256]; best = (long double) astatnum[set/256] * (long double) astatnum[set%256]; bestset = set; deflect = set; for (i=0; i<NUM_PERMS; i++) { reflected = perm(set,i); if (set==reflected) mult++; count = (long double) astatnum[reflected/256] * (long double) astatnum[reflected % 256]; // printf("%ld,%ld->%ld*%ld ", reflected/256, reflected%256, astatnum[reflected/256], astatnum[reflected % 256]); if (count < best) { best = count; a = astatnum[reflected/256]; b = astatnum[reflected%256]; deflect = reflected;} if (reflected < bestset) bestset = reflected; }


// printf("%lo:%d ", set, mult); if (set == bestset) { cost += best; num_apairs++; cost2 += (long double) a * (long double) b; printf(" %ld %ld %ld*%ld = %lu\n", deflect/256, deflect%256, a, b, a*b); } // if (set % 256 == 0) printf("."); } 

  printf("\nTotal cost: %lf\n", (double) cost);
  printf("\nTotal cost alt: %lf\n", (double) cost2);
  printf("\nNumber of a-sets to check: %lf\n", (double) num_apairs);
  printf("\nNumber of Level 1/Level 3 pairs per a-set to check: %lf\n", (double) (cost/num_apairs) );

// init_data();

}

Retrieved from "https://michaelnielsen.org/polymath/index.php?title=Optimal_a-set_pair_code&oldid=1635"