/* Routine to solve Verners 4 cube problem Niels Lund, november 2016 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

int main(int argc, char* argv[])
{
   int i, j, k, index, axnum=0, asum, i1, i2, i3, i0, asum0, asum1, asum2;
   int axn, sol, sidesum0, sidesum1, sidesum2, sidesum3, axphase0, axphase1, axphase2, axphase3;
   int sides00, sides01, sides02, sides03, sides10, sides11, sides12, sides13, sides20, sides21, sides22, sides23; 

   /* blue: 1, green: 2, red: 4, white: 8 */
   int cube[4][6] = {1,2,4,4,8,8, 1,2,2,8,8,4, 1,2,8,4,2,1, 1,2,4,8,4,4};
   int order[3][2][4] = {1,2,4,5, 1,5,4,2, 3,2,0,5, 3,5,0,2, 1,3,4,0, 1,0,4,3};
   int sums[4][3], axes[91][4];

   char sidecolor[4][4][8], color[9][8];

   for (i=0; i<9; i++) sprintf(color[i], "%s", "     ");
   sprintf(color[1], "%s", " blue");
   sprintf(color[2], "%s", "green");
   sprintf(color[4], "%s", "  red");
   sprintf(color[8], "%s", "white");

   for (i=0; i<3; i++) for (j=0; j<4; j++) sums[j][i] = 0;

   printf("\nFind sidesums for each cube for each of three possible axes:\n");
   printf("Colorvalues used;  blue: 1, green: 2, red: 4, white: 8\n");
   k = 0;
   for (j=0; j<4; j++) {	/* For each of the four cubes */
      for (i=0; i<3; i++) {	/* Try three axis directions */
         index = order[i][0][0];
         sums[j][i] += cube[j][index];	/* Add contributions from four sides */
         index = order[i][0][1];
         sums[j][i] += cube[j][index];
         index = order[i][0][2];
         sums[j][i] += cube[j][index];
         index = order[i][0][3];
         sums[j][i] += cube[j][index];
         printf("%2d, cube: %1d, axis: %1d, sidesum: %2d\n", k, j, i, sums[j][i]);
         k++;
      }
   }

   /* we now have know the sums of side face values for each of the 4 cubes and each of three axis orientations */

   printf("\nFind axis combinations for the four cubes where sum of all sidesums is 60:\n");
   axnum = 0;
   k = 0;
   for (i0=0; i0<3; i0++) {
      asum = sums[0][i0];  
      asum0 = asum;
      for (i1=0; i1<3; i1++) {
         asum += sums[1][i1];  
         asum1 = asum;
         for (i2=0; i2<3; i2++) {
            asum += sums[2][i2];  
            asum2 = asum;
            for (i3=0; i3<3; i3++) {
               asum += sums[3][i3];  
               if (asum == 60) {
                  axes[axnum][0] = i0;
                  axes[axnum][1] = i1;
                  axes[axnum][2] = i2;
                  axes[axnum][3] = i3;
                  axnum++;
                  printf("%2d, num: %2d, asum: %2d, axes: %1d %1d %1d %1d\n", k, axnum, asum, i0, i1, i2, i3);
               }
               k++;
               asum = asum2;
            }
            asum = asum1;
         }
         asum = asum0;
      }
   }

   /* we now know which combinations of axis directions for the four cubes which possibly may provide the solution
   int cube[4][6] = {1,2,4,4,8,8, 1,2,2,8,8,4, 1,2,8,4,2,1, 1,2,4,8,4,4};
   int order[3][2][4] = {1,2,4,5, 1,5,4,2, 3,2,0,5, 3,5,0,2, 1,3,4,0, 1,0,4,3}; */

   printf("\nSearch solutions by going through all cube orientations for the %1d valid axis combinations\n", axnum);
   k = 0;
   sol = 0;
   for (axn=0; axn<axnum; axn++) {	/* loop over possible axis combinations */
      for (axphase0=0; axphase0<1; axphase0++) {	/* loop over +/- axis orientation for cube 0 */ 
         i0 = 0;			 /* we do not loop over rotations of cube 0, it can stay fixed */
         index = order[axes[axn][0]][axphase0][(i0+0)%4];
         sidesum0 = cube[0][index];
         sprintf(sidecolor[0][0], "%s", color[cube[0][index]]);
         index = order[axes[axn][0]][axphase0][(i0+1)%4];
         sidesum1 = cube[0][index];
         sprintf(sidecolor[0][1], "%s", color[cube[0][index]]);
         index = order[axes[axn][0]][axphase0][(i0+2)%4];
         sidesum2 = cube[0][index];
         sprintf(sidecolor[0][2], "%s", color[cube[0][index]]);
         index = order[axes[axn][0]][axphase0][(i0+3)%4];
         sidesum3 = cube[0][index];
         sprintf(sidecolor[0][3], "%s", color[cube[0][index]]);
         sides00 = sidesum0;
         sides01 = sidesum1;
         sides02 = sidesum2;
         sides03 = sidesum3;
         for (axphase1=0; axphase1<2; axphase1++) {	/* loop over +/- axis orientation for cube 1 */ 
            sidesum0 = sides00;
            sidesum1 = sides01;
            sidesum2 = sides02;
            sidesum3 = sides03;
            for (i1=0; i1<4; i1++) {			/* loop over rotations of cube 1 */
               index = order[axes[axn][1]][axphase1][(i1+0)%4];
               sidesum0 += cube[1][index];
               sprintf(sidecolor[1][0], "%s", color[cube[1][index]]);
               index = order[axes[axn][1]][axphase1][(i1+1)%4];
               sidesum1 += cube[1][index];
               sprintf(sidecolor[1][1], "%s", color[cube[1][index]]);
               index = order[axes[axn][1]][axphase1][(i1+2)%4];
               sidesum2 += cube[1][index];
               sprintf(sidecolor[1][2], "%s", color[cube[1][index]]);
               index = order[axes[axn][1]][axphase1][(i1+3)%4];
               sidesum3 += cube[1][index];
               sprintf(sidecolor[1][3], "%s", color[cube[1][index]]);
               sides10 = sidesum0;
               sides11 = sidesum1;
               sides12 = sidesum2;
               sides13 = sidesum3;
               for (axphase2=0; axphase2<2; axphase2++) {	/* loop over +/- axis orientation for cube 2 */ 
                  sidesum0 = sides10;
                  sidesum1 = sides11;
                  sidesum2 = sides12;
                  sidesum3 = sides13;
                  for (i2=0; i2<4; i2++) {			/* loop over rotations of cube 2 */
                     index = order[axes[axn][2]][axphase2][(i2+0)%4];
                     sidesum0 += cube[2][index];
                     sprintf(sidecolor[2][0], "%s", color[cube[2][index]]);
                     index = order[axes[axn][2]][axphase2][(i2+1)%4];
                     sidesum1 += cube[2][index];
                     sprintf(sidecolor[2][1], "%s", color[cube[2][index]]);
                     index = order[axes[axn][2]][axphase2][(i2+2)%4];
                     sidesum2 += cube[2][index];
                     sprintf(sidecolor[2][2], "%s", color[cube[2][index]]);
                     index = order[axes[axn][2]][axphase2][(i2+3)%4];
                     sidesum3 += cube[2][index];
                     sprintf(sidecolor[2][3], "%s", color[cube[2][index]]);
                     sides20 = sidesum0;
                     sides21 = sidesum1;
                     sides22 = sidesum2;
                     sides23 = sidesum3;
                     for (axphase3=0; axphase3<2; axphase3++) {	/* loop over +/- axis orientation for cube 3 */ 
                        sidesum0 = sides20;
                        sidesum1 = sides21;
                        sidesum2 = sides22;
                        sidesum3 = sides23;
                        for (i3=0; i3<4; i3++) {			/* loop over rotations of cube 3 */
                           index = order[axes[axn][3]][axphase3][(i3+0)%4];
                           sidesum0 += cube[3][index];
                           sprintf(sidecolor[3][0], "%s", color[cube[3][index]]);
                           index = order[axes[axn][3]][axphase3][(i3+1)%4];
                           sidesum1 += cube[3][index];
                           sprintf(sidecolor[3][1], "%s", color[cube[3][index]]);
                           index = order[axes[axn][3]][axphase3][(i3+2)%4];
                           sidesum2 += cube[3][index];
                           sprintf(sidecolor[3][2], "%s", color[cube[3][index]]);
                           index = order[axes[axn][3]][axphase3][(i3+3)%4];
                           sidesum3 += cube[3][index];
                           sprintf(sidecolor[3][3], "%s", color[cube[3][index]]);
                           if ((sidesum0 == 15) && (sidesum1 == 15) && (sidesum2 == 15) && (sidesum3 == 15)) { 
                              sol++;
printf("Solution %2d: cube0: ax %2d rot %1d, cube1:  ax %2d rot %1d, cube2: ax %2d rot %1d,cube3:  ax %2d rot %1d\n",
                              sol, (1-2*axphase0)*axes[axn][0], i0, (1-2*axphase1)*axes[axn][1], i1, 
                              (1-2*axphase2)*axes[axn][2], i2, (1-2*axphase3)*axes[axn][3], i3);
                              printf("         cube0 cube1 cube2 cube3\n");
                              printf(" side 0: %5s %5s %5s %5s\n", 
                                 sidecolor[0][0], sidecolor[1][0], sidecolor[2][0], sidecolor[3][0]);
                              printf(" side 1: %5s %5s %5s %5s\n", 
                                 sidecolor[0][1], sidecolor[1][1], sidecolor[2][1], sidecolor[3][1]);
                              printf(" side 2: %5s %5s %5s %5s\n", 
                                 sidecolor[0][2], sidecolor[1][2], sidecolor[2][2], sidecolor[3][2]);
                              printf(" side 3: %5s %5s %5s %5s\n", 
                                 sidecolor[0][3], sidecolor[1][3], sidecolor[2][3], sidecolor[3][3]);
                           } 
                           k++;
                           sidesum0 = sides20;
                           sidesum1 = sides21;
                           sidesum2 = sides22;
                           sidesum3 = sides23;
                        }
                     }
                     sidesum0 = sides10;
                     sidesum1 = sides11;
                     sidesum2 = sides12;
                     sidesum3 = sides13;
                  }
               }
               sidesum0 = sides00;
               sidesum1 = sides01;
               sidesum2 = sides02;
               sidesum3 = sides03;
            }
         }
      }
   }
   printf("Finished. Found %1d solution, tried %5d cube orientations\n\n", sol, k);
}