//check up program to prove no 12x12 franklin magic square exists
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define abs(X) ((X)>0?(X):(-(X)))
int nocand[7];
int cand[7][144];
char candinv[7][144];
int x[12],y[12];;
char sel[144];
int sq[12][12];
int wcnt=0,acnt=0,fcnt=0;

int unsel(int j) { sel[j]=0; }

int cleanrow(int rwnr, int numf) {
  int i,nf2;
  nf2=2*numf;
  for (i=0;i<nf2;i+=2) sel[ sq[rwnr][i] ]=0;
}

int cleanfullrow(int rwnr, int numf) {
  int i;
  for (i=0;i<numf;i++) sel[ sq[rwnr][i] ]=0;
}

int cleancol(int clnr, int numf) {
  int i;
  // only for odd cols
  for (i=0;i<numf && i<6 ;i++) sel[ sq[2*i][clnr] ]=0;
  for (i=6;i<numf ;i++) sel[ sq[2*i-11][clnr] ]=0;
}

int addcol(int clnr, int ffld) {
  int i,k,kk;
  // only for odd cols
  for (i=0;i<12;i+=2) {
   k = y[i]+ffld;
   kk = 143-k;
   if (k<144 && k>=0 && sel[kk]==0) {
    sel[kk]=1;
    sq[i][clnr] = kk;
   }
   else {
    break;
   }
  }
  if (i<12) return (i>>1);
  for (i=6;i<9;i++) {
   k = y[2*i-11]+ffld;
   if (k<144 && k>=0 && sel[k]==0) {
    sel[k]=1;
    sq[2*i-11][clnr] = k;
   }
   else {
    break;
   }
  }
  return(i);
}

int addfullrow(int rwnr, int ffld) {
  int i,k,kk;
  for (i=0;i<12;i++) {
   k = x[i]+ffld;
   kk = ( (rwnr+i)%2==0 ? k : 143-k );
   if (k<144 && k>=0 && sel[kk]==0) {
    sel[kk]=1;
    sq[rwnr][i] = kk;
   }
   else {
    break;
   }
  }
  return (i);
}

int addevenrow(int rwnr, int ffld) {
  int i,k;
  // for even rwnr only
  for (i=0;i<12;i+=2) {
   k = x[i]+ffld;
   if (k<144 && sel[k]==0) {
    sel[k]=1;
    sq[rwnr][i] = k;
   }
   else {
    break;
   }
  }
  return( (i>>1) );
}

int addoddrow(int rwnr, int ffld) {
  int i,k;
  // for odd rwnr only
  for (i=0;i<12;i+=2) {
   k = 143-x[i]-ffld;
   if (k>=0 && sel[k]==0) {
    sel[k]=1;
    sq[rwnr][i] = k;
   }
   else {
    break;
   }
  }
  return( (i>>1) );
}

int extend6x6();

int main(int argc, char *argv[]) {
  int bmax,Tmax,Tcomp,x2,x4,x6,x8,x10,y2,y4,y6,y8,y10,f2,f4,f6,f8,f10,\
      lb6,ub6,lb8,ub8,lby2,y2cnt,y4cnt,y6cnt,y8cnt,y10cnt,\
      lby6,lby8,uby8,lby10,uby10,i,j,k,T;
  int lb2,ub2,res2;
  double cnt66=0.0, cand66=0.0;
  int xcnt = 0,ycnt=0,zcnt=0,vcnt=0;

  memset(sel,0,144);
  for (i=0;i<12;i++) for (j=0;j<12;j++) sq[i][j]=0;

  lb2 = 1; ub2 = 70;
  if (argc>1) {res2=atoi(argv[1]); lb2=ub2=res2; }
  x[0]=0; y[0]=0;
  sel[0]=1;
  sq[0][0]=0;
  bmax = 71;
  for (x2=lb2; x2<=ub2; x2++) {
   x[2]=x2;
   sel[x2]=1;
   sq[0][2]=x2;
   for (x4=x2+1; x4<=bmax; x4++) {
    x[4]=x4;
    sel[x4]=1;
    sq[0][4]=x4;
    ub6 = (x2+x4-3)/3;
    for (x6=1; x6<=ub6; x6++) if (x6!=x2 && x6!=x4) {
     x[6]=x6;
     sel[x6]=1;
     sq[0][6]=x6;
     lb8 = (x6+1>x2+x4-x6-bmax?x6+1:x2+x4-x6-bmax);
     ub8 = (x2+x4-x6-1)>>1;
     for (x8=lb8; x8<=ub8; x8++) if (x8!=x2 && x8!=x4) {
      x[8]=x8;
      sel[x8]=1;
      sq[0][8]=x8;
      x10 = x2+x4-x6-x8;
      if (x10!=x2 && x10!=x4) {
       x[10]=x10;
       sel[x10]=1;
       sq[0][10]=x10;
xcnt++; 
       Tmax = (x4>x10?x4:x10);
       Tcomp = 143 - Tmax;

       memset(candinv[0],0,144);
       memset(candinv[0],1,Tcomp+1);
       for (i=0;i<12;i+=2) candinv[0][x[i]]=0;
       for (i=2;i<10;i+=2) for (j=i+2;j<12;j+=2)
        candinv[0][abs(x[j]-x[i])]=0;
       for (i=0,j=0;i<=Tcomp;i++) if (candinv[0][i]) cand[0][j++]=i;
       nocand[0]=j;
       // filling 5 more rows distinguish between y4<>Tmax
       for (y4cnt=0; y4cnt<nocand[0]; y4cnt++) {
        y4 = cand[0][y4cnt];
        if (y4 >= ((Tmax+6)>>1) ) break;
       }
       for ( ; y4cnt<nocand[0]; y4cnt++) {
        y4 = cand[0][y4cnt];
        if (y4>Tcomp) {
         y4cnt=nocand[0]; continue;
        }
        y[4]=y4;
        f4 = addevenrow(4,y4);
        if (f4<6) {
         cleanrow(4,f4);
         continue;
        }
ycnt++;
        lby2 = (y4<Tmax?Tmax+4-y4:1);
        for (y2cnt=0; y2cnt<nocand[0]; y2cnt++) {
         y2 = cand[0][y2cnt];
         if (y2>=lby2) break;
        }
        for (y2cnt=0; y2cnt<nocand[0]; y2cnt++) {
         y2 = cand[0][y2cnt];
         if (y2>=y4) {
          y2cnt = nocand[0]; continue;
         }
         y[2]=y2;
         f2 = addevenrow(2,y2);
         if (f2<6) {
          cleanrow(2,f2);
          continue;
         }
zcnt++; 
         memset(candinv[1],0,144);
         for (i=0,j=0;i<nocand[0];i++) {
          T = cand[0][i];
          for (k=0;k<12;k+=2) {
           if (sel[T+x[k]]==1) break;
          }
          if (k>=12) {
           candinv[1][T]=1;
           cand[1][j++]=T;
          }
         }
         nocand[1]=j;

         uby10 = (y2+y4-3<Tcomp?y2+y4-3:Tcomp);
         lby10 = (y2+y4+2)/3;
         if (y4<Tmax && lby10<=Tmax) lby10=Tmax+1;
         for (y10cnt=nocand[1]-1; y10cnt>=0 ; y10cnt--) {
          y10 = cand[1][y10cnt];
          if (y10<=uby10) break;
         }
         for ( ; y10cnt>1 ; y10cnt--) {
          y10 = cand[1][y10cnt];
          if (y10<lby10) {
             y10cnt=0; continue;
          }
          y[10]=y10;
          f10 = addevenrow(10,y10);
          if (f10<6) {
           cleanrow(10,f10);
           continue;
          }
vcnt++;
          uby8 = y2+y4-y10-1;
          lby8 = (y2+y4-y10+1)>>1;
          for (y8cnt=y10cnt-1; y8cnt>0; y8cnt--) {
           y8 = cand[1][y8cnt];
           if (y8<=uby8) break;
          }
          for ( ; y8cnt>0; y8cnt--) {
           y8 = cand[1][y8cnt];
           if (y8<lby8) {
              y8cnt=0; continue;
           }
           y[8]=y8;
           f8 = addevenrow(8,y8);
           if (f8<6) {
            cleanrow(8,f8);
            continue;
           }
           y6 = y2+y4-y10-y8;
           if (candinv[1][y6]==1) {
            y[6]=y6;
            f6 = addevenrow(6,y6);
            if (f6<6) {
             cleanrow(6,f6);
             cleanrow(8,f8);
             continue;
            }
            // now 6x6 subsquare has been filled, proceed by
            cand66 += 1.0;
                         ////////////
            cnt66 += 0.0+extend6x6();
                         ////////////
            cleanrow(6,6);
           } //end y6
           cleanrow(8,6);
          } //end y8
          cleanrow(10,6);
         } //end y10
         cleanrow(2,6);
        } //end y2 
        cleanrow(4,6);
       } //end y4
       unsel(x10);
      } //end if x10
      unsel(x8);
     } //end x8
     unsel(x6);
    } //end x6
    unsel(x4);
   } //end x4
   unsel(x2);
  } //end x2
  printf("counting %d tuples x2..10, and\n %d, %d, %d fitting y4, y2, y10\n",
         xcnt,ycnt,zcnt,vcnt);
  printf("counted %10.0lf 6 by 6 sols out of\n %10.0lf cands for x2 = %d\n",
         cnt66, cand66, res2);
  printf("found %d nine by nine-s, %d of which with potential\n",wcnt,acnt);
  printf("detected %d Franklin Magic Squares\n", fcnt);
}

int tripcnt(int nrc, int *cnd, char *invcnd, int sum) {
  // check if nr of triples a<b<c in cnd with a+b+c=sum EXCEEDS 10
  int i,j,k,cnt=0;
  for (i=0;i<nrc-2 && cnd[i]<=(sum-3)/3 ;i++)
  for (j=i+1;j<nrc-1 && cnd[j]<=( (sum-cnd[i]-1)>>1) ;j++) {
   k = sum-cnd[i]-cnd[j];
   if (k<144 && invcnd[k]==1 && ++cnt>=10) return (cnt);
  }
  return (cnt);
}

int pas1=0,pas2=0,pas3=0,pas4=0,pas5=0;

int extend9x9();

int extend6x6() {
  int i,j,k,m,f,Ymax,Ycomp,Xmax,Xcomp;
  int xf1,xf3,xf5,f1,f3,f5,uby1,ubx1,T;
  int x1cnt,x3cnt,y1cnt,y3cnt;
  int x1,x3,x5,y1,y3,y5,ubx3,uby3;

  // compute odd column leader candidates
  Ymax=(y[4]>y[10]?y[4]:y[10]);
  Ycomp = 143-Ymax;
  memset(candinv[2],0,144);
  for (i=0,j=0;i<=Ycomp;i++) {
   for (k=0; k<12; k+=2) {
    if (sel[143-i-y[k]]==1) break;
   }
   if (k>=12) {
    candinv[2][i] = 1;
    cand[2][j++]=i;
   }
  }
  nocand[2]=j;
  if (nocand[2]<6) return (0);
  if (tripcnt(nocand[2],cand[2],candinv[2],x[2]+x[4]) < 2) return(0);
  // compute odd row leader candidates
  Xmax=(x[4]>x[10]?x[4]:x[10]);
  Xcomp = 143-Xmax;
  memset(candinv[3],0,144);
  for (i=0,j=0;i<=Xcomp;i++) {
   for (k=0; k<12; k+=2) {
    if (sel[143-i-x[k]]==1) break;
   }
   if (k>=12) {
    candinv[3][i] = 1;
    cand[3][j++]=i;
   }
  }
  nocand[3]=j;
  if (nocand[3]<6) return (0);
  if ( (f=tripcnt(nocand[3],cand[3],candinv[3],y[2]+y[4])) < 2){
     return(0);
  }
  
  uby1=(y[2]+y[4]-3)/3;
  for (y1cnt=0; y1cnt<nocand[3]-2; y1cnt++) {
   y1 = cand[3][y1cnt];
   if (y1>uby1) {
    y1cnt = nocand[3];
    continue;
   }
   y[1]=y1;
   f1 = addoddrow(1,y1);
   if (f1<6) {
    cleanrow(1,f1);
    continue;
   }
   uby3=(y[2]+y[4]-y1-1)>>1;
   for (y3cnt=y1cnt+1; y3cnt<nocand[3]-1; y3cnt++) {
    y3 = cand[3][y3cnt];
    if (y3>uby3) {
     y3cnt = nocand[3];
     continue;
    }
    y[3]=y3;
    f3 = addoddrow(3,y3);
    if (f3<6) {
     cleanrow(3,f3);
     continue;
    }
    y5 = y[2]+y[4]-y1-y3;
    if (y5<144 && candinv[3][y5]==1) {
     y[5]=y5;
     f5 = addoddrow(5,y5);
     if (f5<6) {
      cleanrow(5,f5);
      cleanrow(3,f3);
      continue;
     }
 
//first sharpen candidate list for col leaders: cand2 -> cand4
     //for (i=0;i<144;i++) candinv[4][i]=candinv[2][i];
     memset(candinv[4],0,144);
     for (i=0,j=0;i<nocand[2];i++) {
      T = cand[2][i];
      for (m=0;m<12;m+=2) {
       if (sel[143-T-y[m]]==1) break;
      }
      if (m>=12) {
       for (k=1;k<7;k+=2) {
        if(T+y[k]>143 || sel[T+y[k]]==1) break;
       }
       if (k>=7) {
        candinv[4][T]=1;
        cand[4][j++]=T;
       }
      }
     }
     nocand[4]=j;
     if ( nocand[4]<6 ||
          tripcnt(nocand[4],cand[4],candinv[4],x[2]+x[4]) < 2 ){
      cleanrow(5,f5);
      cleanrow(3,f3);
      continue;
     }

     ubx1=(x[2]+x[4]-3)/3;
     for (x1cnt=0; x1cnt<nocand[4]-2; x1cnt++) {
      x1 = cand[4][x1cnt];
      if (x1>ubx1) {
       x1cnt = nocand[4];
       continue;
      }
      x[1]=x1;
      xf1 = addcol(1,x1);
      if (xf1<9) {
       cleancol(1,xf1);
       continue;
      }
      ubx3=(x[2]+x[4]-x1-1)>>1;
      for (x3cnt=x1cnt+1; x3cnt<nocand[4]-1; x3cnt++) {
       x3 = cand[4][x3cnt];
       if (x3>ubx3) {
        x3cnt = nocand[4];
        continue;
       }
       x[3]=x3;
       xf3 = addcol(3,x3);
       if (xf3<9) {
        cleancol(3,xf3);
        continue;
       }
       x5 = x[2]+x[4]-x1-x3;
       if (candinv[4][x5]==1) {
        x[5]=x5;
        xf5 = addcol(5,x5);
        if (xf5<9) {
         cleancol(5,xf5);
         cleancol(3,xf3);
         continue;
        }

wcnt += 1; 
        ////////////
        extend9x9();
        ////////////
        cleancol(5,9);
       } //end x5
       cleancol(3,9);
      } //end x3
      cleancol(1,9);
     } //end x1
     cleanrow(5,6);
    } //end y5
    cleanrow(3,6);
   } //end y3
   cleanrow(1,6);
  } //end y1
  return (1);
  
}

int extend9x9() {
  int i,j,k,m,T, x7cnt,x9cnt,y7cnt,y9cnt;
  int x7,x9,x11,y7,y9,y11,ubx7,ubx9,uby7,uby9,xf7,xf9,xf11,yf7,yf9,yf11;
//first sharpen candidate list for col leaders: cand4 -> cand6
  memset(candinv[6],0,144);
  for (i=0,j=0;i<nocand[4];i++) {
   T = cand[4][i];
   for (m=0;m<12;m+=2) {
    if (sel[143-T-y[m]]==1) break;
   }
   if (m>=12) {
    for (k=1;k<7;k+=2) {
     if(sel[T+y[k]]==1) break;
    }
    if (k>=7) {
     candinv[6][T]=1;
     cand[6][j++]=T;
    }
   }
  }
  nocand[6]=j;
  if ( nocand[6]<3 ||
       tripcnt(nocand[6],cand[6],candinv[6],x[2]+x[4]) < 1 ){
   return (0);
  }

//next sharpen candidate list for row leaders: cand3 -> cand5
  memset(candinv[5],0,144);
  for (i=0,j=0;i<nocand[3];i++) {
   T = cand[3][i];
   for (m=0;m<12;m+=2) {
    if (sel[143-T-x[m]]==1) break;
   }
   if (m>=12) {
    for (k=1;k<7;k+=2) {
     if(T+x[k]>143 || sel[T+x[k]]==1) break;
    }
    if (k>=7) {
     candinv[5][T]=1;
     cand[5][j++]=T;
    }
   }
  }
  nocand[5]=j;
  if ( nocand[5]<3 ||
       tripcnt(nocand[5],cand[5],candinv[5],y[2]+y[4]) < 1 ){
   return (0);
  }

  acnt++;
  if (acnt==1) {
   for (i=0;i<12;i++) {
    for (j=0;j<12;j++) printf(" %3d", sq[i][j]);
    printf("\n");
   }
   fflush(stdout);
  }
  // now trying for real

  ubx7=(x[2]+x[4]-3)/3;
  for (x7cnt=0; x7cnt<nocand[6]-2; x7cnt++) {
   x7 = cand[6][x7cnt];
   if (x7>ubx7) {
    x7cnt = nocand[6];
    continue;
   }
   x[7]=x7;
   xf7 = addcol(7,x7);
   if (xf7<9) {
    cleancol(7,xf7);
    continue;
   }
   ubx9=(x[2]+x[4]-x7-1)>>1;
   for (x9cnt=x7cnt+1; x9cnt<nocand[6]-1; x9cnt++) {
    x9 = cand[6][x9cnt];
    if (x9>ubx9) {
     x9cnt = nocand[6];
     continue;
    }
    x[9]=x9;
    xf9 = addcol(9,x9);
    if (xf9<9) {
     cleancol(9,xf9);
     continue;
    }
    x11 = x[2]+x[4]-x7-x9;
    if (candinv[6][x11]==1) {
     x[11]=x11;
     xf11 = addcol(11,x11);
     if (xf11<9) {
      cleancol(11,xf11);
      cleancol(9,xf9);
      continue;
     }
// finally fill last three rows 7,9,11:
     uby7=(y[2]+y[4]-3)/3;
     for (y7cnt=0; y7cnt<nocand[5]-2; y7cnt++) {
      y7 = cand[5][y7cnt];
      if (y7>uby7) {
       y7cnt = nocand[5];
       continue;
      }
      y[7]=y7;
      yf7 = addfullrow(7,y7);
      if (yf7<12) {
       cleanfullrow(7,yf7);
       continue;
      }
      uby9=(y[2]+y[4]-y7-1)>>1;
      for (y9cnt=y7cnt+1; y9cnt<nocand[5]-1; y9cnt++) {
       y9 = cand[5][y9cnt];
       if (y9>uby9) {
        y9cnt = nocand[5];
        continue;
       }
       y[9]=y9;
       yf9 = addfullrow(9,y9);
       if (yf9<12) {
        cleanfullrow(9,yf9);
        continue;
       }
       y11 = y[2]+y[4]-y7-y9;
       if (y11<144 && candinv[5][y11]==1) {
        y[11]=y11;
        yf11 = addfullrow(11,y11);
        if (yf11<12) {
         cleanfullrow(11,yf11);
         cleanfullrow(9,yf9);
         continue;
        }
  
        // if here found FMS
        fcnt++;

        printf("Found Franklin Magic 12 by 12 Square\n");
        for (i=0; i<12;i++) {
         for (j=0;j<12;j++) printf(" %3d",sq[i][j]);
         printf("\n");
        }

        cleanfullrow(11,12);
       } //end y11
       cleanfullrow(9,12);
      } //end y9
      cleanfullrow(7,12);
     } //end y7
     cleancol(11,9);
    } //end x11
    cleancol(9,9);
   } //end x9
   cleancol(7,9);
  } //end x7

  return (1);
}