The number of magic squares of order 6

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
<file c ms-20230905-m.c>
/*
 
  See https://magicsquare6.net/ .
 
  Updated on 2023/09/05.
  Authored by Hidetoshi Mino hidetoshimino@gmail.com .
 
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <pthread.h>
#ifndef noMD5
  #include <openssl/md5.h>
#endif
#define MAX( x, y ) ( (x) > (y) ? (x) : (y) )

#ifndef N
  #define N 6
#endif

#ifndef NTH
  #define NTH 1
#endif

#define NN (N*N)
#define AllBits ( ( 1ULL << NN ) - 1 )

typedef long long unsigned u64;

inline int topbit( u64 v ) {

   v = v & ( ~( v >> 1 ) );
   float f32 = (float) v;
   int pos = *( (int *) &f32 );
   pos =( pos >> 23 ) - 0x7f;
   return pos;
/*
  int pos = 0;
  int v32 = v >> 32;
  if( v32 ) {
      pos = 32;
  } else {
     v32 = (int) v;
  }
  if( v32 & 0xffff0000UL ) { pos +=16;  v32 >>= 16; }
  if( v32 & 0x0000ff00UL ) { pos += 8;  v32 >>=  8; }
  if( v32 & 0x000000f0UL ) { pos += 4;  v32 >>=  4; }
  if( v32 & 0x0000000cUL ) { pos += 2;  v32 >>=  2; }
  if( v32 & 0000000002UL ) pos++;
  return pos;
*/
}

u64 reverse64( u64 x ) {
   u64 y;
   y = ( ( x & 0x5555555555555555ULL ) << 1 );
   x = ( y | ( ( x >> 1 ) & 0x5555555555555555ULL ) );
   y = ( ( x & 0x3333333333333333ULL ) << 2 );
   x = ( y | ( ( x >> 2 ) & 0x3333333333333333ULL ) );
   y = ( ( x & 0x0F0F0F0F0F0F0F0FULL ) << 4 );
   x = ( y | ( ( x >> 4 ) & 0x0F0F0F0F0F0F0F0FULL ) );
   y = ( ( x & 0x00FF00FF00FF00FFULL ) << 8 );
   x = ( y | ( ( x >> 8 ) & 0x00FF00FF00FF00FFULL ) );
   y = ( ( x & 0x0000FFFF0000FFFFULL ) << 16 );
   x = ( y | ( ( x >> 16 ) & 0x0000FFFF0000FFFFULL ) );
   y = ( ( x & 0x00000000FFFFFFFFULL ) << 32 );
   x = ( y | ( ( x >> 32 ) & 0x00000000FFFFFFFFULL ) );
   return x;
}

int check_topbit( ) {

   for( int i = 0; i < NN; i ++ ) {
      // printf( "%d %09llx %d\n", i, 1ULL << i, topbit( 1ULL << i ) );
      if( topbit( 1ULL << i ) != i ) return 1;
   }
   for( int i = 1; i < NN; i ++ ) {
      // printf( "%d %09llx %d\n"
      //         , i, ( 1ULL << i ) - 1, topbit( ( 1ULL << i ) - 1 ) );
      if( topbit( ( 1ULL << i ) - 1 ) != i - 1 ) return 1;
   }
   fprintf( stderr, "topbit() works correctly.\n" );
   return 0;
}

int factorial( int i ) {
   if( i < 1 ) return 1;
   return factorial( i - 1 ) * i;
}

#define OW 3   // octal bit witdh
#define OMASK ( ( 1 << OW ) - 1 )  // octal bit mask
      
int pcode( int index, int n )
{  
  int posv = 076543210; // octal number
  int onev = 011111111; // octal number
  int res = 0; 
  n--;
  for(  ; n > 0; n-- ) {
     int pos = ( index >> ( n * OW ) ) & OMASK ;
     int shift = pos * OW;
     pos = ( ( posv >> shift ) & OMASK );
     shift += OW;
     posv -= ( ( onev >> shift ) << shift );
     res += pos;
     res *= n;
  }
  return res;
}

int pinv( int pcode, int n )
{
   int res = 0;
   for( int b = 1; b < n; b ++ ) {
      int p = pcode % ( b + 1 );
      int shift = p * OW;
      int upper = ( ( res >> shift ) << ( shift + OW ) );
      res = upper | ( b << shift ) | ( res & ( ( 1 << shift ) - 1 ) );
      pcode /= ( b + 1 );
   }
   return res;
}

int inv_indx ( int index )
{
  int res = 0;
  int pos = 0;
  for(  ; index != 0; index >>= OW ) {
     res |= ( pos << ( ( index & OMASK ) * OW ) );
     pos ++;
  }

  return res;
}

void GenPerm( int order, const int mask_w, int* res, int res_size )
{
   if( order < 2 ) {
      *res = 0;
      return;
   }

   int sub_size = res_size / order;
   int sub_res[ sub_size ];

   GenPerm( order - 1, mask_w, sub_res, sub_size  );

   unsigned up_mask  = ( 1u << ( order - 1 ) * mask_w ) - 1;
   unsigned low_mask = up_mask;
   unsigned insert = ( order - 1 ) <<  ( order - 1 ) * mask_w;
   for( int i = 0; i < order; i ++ ) {
      for( int j = 0; j < sub_size; j ++ ) { 
         int upper = sub_res[ j ] & ( up_mask ^ low_mask );
         int lower = sub_res[ j ] & low_mask;
         *(res++) = ( upper << mask_w ) | insert | lower ;
      }
      low_mask >>= mask_w;
      insert  >>= mask_w;
   }
}

/* = globals ======================= */

#define MPERM 720
#define MPERMW ( ( MPERM + 31 ) / 32 )
unsigned diag_chk_a[ MPERM ][ MPERMW ];

#define MVECS 32768
int nvecs = 0;
u64 vecs[ MVECS ];
u64 vec_elem[ MVECS ]; // Element numbers are scrumbled by (2^n) % 37.

u64 count_g;

/*= Thread shared  ========================= */

// read/write

int* third_row;
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;

// read only

int max_row, second_row;
u64 rev_max;  // Reverse( vecs[ max_row ] )

int*    rows;  int nrows;
u64* cols;  int ncols;
u64* diags; int ndiags;

/* = Thread local ========================== */

__thread u64 col_vecs[ N ]; 
__thread u64 row_vecs[ N ];
__thread int which_row[ 37 ]; 
                         // 37 is the common magic number for order <= 6.
__thread unsigned sub_count;

__thread  int*    vr;
__thread  u64* vc;
__thread  u64* vx;

/* ========================================= */
void Count( );
void MakeRows( int level, u64 rem, int* vrb, int* vre,
        u64* vcb, u64* vce, u64* vxb, u64* vxe, int single, u64* countp );
void ThirdRow( );
void MakeCols( int level, u64 rem, 
        u64* vcb, u64* vce, u64* vxb, u64* vxe, int single );
void MakeColsLast( int level, u64 rem, 
        u64* cols,          u64* vxb, u64* vxe, int single );
unsigned MakeDiag( unsigned* diags, int ndiags );
int DiagPos( u64* vxb, u64* vxe, u64 col_vec, unsigned* cxe );

int Cross( int v1, int v2 ) {
   u64 v = vecs[v1] & vecs[v2];
   return v && ( ! ( v & ( v - 1 ) ) );
}

u64 Reverse( u64 v ) { return reverse64( v ) >> (64 - NN); }

void MakeDiagCheck( int n )
{
   int perm_size = factorial( n );
   int perm[ perm_size ];
   static const int mask_w = 3;
   GenPerm( n, mask_w, perm, perm_size );

   unsigned lower_mask = ( 1 << mask_w ) - 1;
   unsigned pos_mask = ( 1 << ( mask_w * (n-1) ) ) - 1;

   for( int dx = 0; dx < perm_size; dx ++ ) {
      unsigned* va = diag_chk_a[ pcode( ( perm[ dx ] & pos_mask ) << 3, N ) ];
      for( int j = 0; j < MPERMW; j ++ ) {
         va[ j ] = 0U;
      }
      int c[ n ];
      for( int k = 0; k < n; k ++ ) {
         c[ ( perm[ dx ] >> ( (n-1-k)*mask_w ) ) & lower_mask ] = k;
      }
      for( int cx = 0; cx < perm_size; cx ++ ) {
         int r[ n ];
         int ncross = 0;
         for( int k = 0; k < n; k ++ ) {
            r[ k ] = ( perm[ cx ] >> ( (n-1-k)*mask_w ) ) & lower_mask;
            if( c [ r[ k ] ] == k ) ncross ++;
         }
         if( ncross == N % 2 &&
             r[ c[ r[ c[ 0 ] ] ] ] == 0 &&
             r[ c[ r[ c[ 1 ] ] ] ] == 1 &&
             r[ c[ r[ c[ 2 ] ] ] ] == 2 )
         {
            int pos_code = pcode( ( perm[ cx ] & pos_mask ) << 3, N );
            va[ pos_code >> 5 ] |= 1U << ( pos_code & ( ( 1 << 5 ) - 1 ) );
         }
      }
   }
}

void GenVecs( int max_ent, int rem_sum, int level, u64 vp ) 
{
  static int elem[ N ];

  if ( level >= 0 ) {
    for ( int ent = max_ent; ent > 0 ; --ent ) {
      u64 v = 1ULL << ( ent - 1 );
      elem[ level ] = v % 37;    // A trick to trasform 2^n into [1:36]
      if ( rem_sum >= ent ) 
        GenVecs( ent - 1, rem_sum - ent, level - 1, vp | v );
    }
  }
  else {
    if ( rem_sum == 0 )  {
       vecs[ nvecs ] = vp;
       u64 ve = 0ULL;
       for( int i = 0; i < N; i ++ ) {
          ve <<= 6;  // 2^6 >= 37
          ve |= elem[ i ];
       }
       vec_elem[ nvecs ] = ve;
       nvecs ++;
    }
  }
  if( nvecs > MVECS ) {
    fprintf( stderr, "nvecs > MVECS\n" );
    exit(1);
  }
}

void Count( ) 
{

  rows  = malloc( sizeof( int ) * nvecs );
  cols  = malloc( sizeof( u64 ) * nvecs );
  diags = malloc( sizeof( u64 ) * nvecs );

  rev_max = Reverse( vecs[ max_row ] );

  int* rowp = rows;
  for( int i = second_row + 1; i < nvecs; ++ i ) 
    if(        ( ! ( vecs[ max_row  ] & vecs[ i ] ) ) 
           && ( ! ( vecs[ second_row ] & vecs[ i ] ) )
           && vecs[ max_row ] >= Reverse( vecs[ i ] ) )  
    {
       *(rowp++) = i;
    }
  nrows = rowp - rows;
  // printf(" %d\n", nrows );

  u64* colp = cols;
  for( int i = max_row + 1; i < nvecs; ++ i ) 
    if( Cross( max_row, i ) && Cross( second_row, i ) 
        && vecs[ max_row ] >= Reverse( vecs[ i ] ) )
    {
       *(colp++) = vecs[ i ];
    }
  ncols = colp - cols;

  u64* diagp = diags;
  for( int i = 0; i < nvecs; ++ i ) 
    if( Cross( max_row, i ) && Cross( second_row, i ) ) 
       *(diagp++) = vecs[ i ];
  ndiags = diagp - diags;

  third_row = rows;

  u64 count[ NTH ];
  pthread_t thread[ NTH ];

  for( int i = 0; i < NTH; i ++ ) {
     count[ i ] = 0ULL;
     if (pthread_create( thread + i , NULL, (void*) &ThirdRow, count + i ) ) {
        fprintf( stderr, "pthread_create error! \n" );
        exit( 1 );
     }
  }

  for( int i = 0; i < NTH; i ++ ) {
     if( pthread_join( thread[ i ], NULL ) ) {
        fprintf( stderr, "pthread_join error! \n" );
        exit( 1 );
     }
     count_g += count[ i ];
  }
  free( diags );
  free( cols );
  free( rows );

  return;
}

void ThirdRow( u64* countp )
{
  row_vecs[ N - 1 ] = vecs[ max_row ];
  row_vecs[ N - 2 ] = vecs[ second_row ];

  // 37 is the common magic number for order <= 6.
  for( int j = 0; j < 37; ++ j ) which_row[ j ] = 0; //magic index range !!
  u64 r1 = vec_elem[ max_row ];
  u64 r2 = vec_elem[ second_row ];
  for( int j = 0; j < N; ++ j ) {
     which_row[ r1 & 0x3f ] = N - 1;
     which_row[ r2 & 0x3f ] = N - 2;
     r1 >>= 6;
     r2 >>= 6;
  }

  int vrsize = nrows  * MAX( N - 4, 1 );
  int vcsize = ncols  * MAX( 2*N - 5, 2 ); 
  int vxsize = ndiags * MAX( 2*N - 4, 3 );
  vr = malloc( sizeof( int ) * vrsize );
  vc = malloc( sizeof( u64 ) * vcsize );
  vx = malloc( sizeof( u64 ) * vxsize );
  vr[ vrsize - 1 ] = 0U;
  vc[ vcsize - 1 ] = 0ULL;
  vx[ vxsize - 1 ] = 0ULL;

  int single = rev_max == vecs[ max_row ] || rev_max == vecs[ second_row ];
  u64 rem = AllBits ^ vecs[ max_row ] ^ vecs[ second_row ]; 
  u64 must_bit = ( 1ULL <<  topbit( rem ) );

  int* vrb;

  for(   ; /* vrb <  rows + nrows - N + 3 */; /* vrb ++ */ ) {

     if( pthread_mutex_lock( &mutex ) ) {
        fprintf( stderr, "Mutex lock error!\n" );
        exit( 1 );
     }
        vrb = third_row;
        third_row ++;

     if( pthread_mutex_unlock( &mutex ) ) {
        fprintf( stderr, "Mutex unlock error!\n" );
        exit( 1 );
     }

     if( vrb > rows + nrows - N + 2 ) break;

     u64 row_vec = vecs[ *vrb ];
     if( ! ( row_vec & must_bit ) ) break;

     u64* vcp = cols;
     u64* vcq = vc;
     while( vcp < cols + ncols ) {
        u64 x = row_vec & ( *vcp );
        *vcq = *vcp;
        vcq += ( ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 ) );
        vcp ++;
     }
     if( vcq - vc < N ) continue;

     u64* vxp = diags;
     u64* vxq = vx;
     while( vxp < diags + ndiags ) {
        u64 x = row_vec & ( *vxp );
        *vxq = *vxp;
        vxq += ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 );
        vxp ++;
     }
     if ( vxq - vx < 2 ) continue;

     row_vecs[ N - 3 ] = row_vec;
     u64 ve = vec_elem[ *vrb ];
     for( int j = 0; j < N; ++ j ) {
        which_row[ ve & 0x3f ] = N - 3;
        ve >>= 6;
     }

     if( N > 4 ) {

        int* vrp = vrb + 1;
        int* vrq = vr;
        while( vrp < rows + nrows ) {
           *vrq = *vrp;
           vrq += ( row_vec & ( vecs[ *vrp ] ) ) ? 0 : 1;
           vrp ++;
        }
        // if( vrq - vre < N - 3 ) continue; // Don't use continue here!!!

#if N > 5
#ifdef V
        fprintf( stderr,
                " level: %d %09llx  %llu \n", 3, row_vec, *countp );
#endif
#endif
        MakeRows( N - 4, rem ^ row_vec, vr, vrq, vc, vcq, vx, vxq, 
                                   single || ( rev_max == row_vec ), countp );
     } else {
#if N > 3
        row_vecs[ 0 ] = rem ^ row_vec;
        if( ( ! ( row_vecs[ 0 ] & 1 ) ) || 
              row_vecs[ N - 1 ] >= Reverse( row_vecs[ 0 ] ) ) 
#endif
        {  
          sub_count = 0U;
            MakeCols( N - 1, AllBits, vc, vcq, vx, vxq, 
                    single || rev_max == row_vec || rev_max == row_vecs[ 0 ] ); 
          *countp += sub_count;
        }
     }
     ve = vec_elem[ *vrb ];
     for( int j = 0; j < N; ++ j ) {
        which_row[ ve & 0x3f ] = 0;
        ve >>= 6;
     }
  }
  if ( vr[ vrsize - 1 ] != 0U || 
       vc[ vcsize - 1 ] != 0ULL ||
       vx[ vxsize - 1 ] != 0ULL ) {
     fprintf( stderr, " Array over run ! \n");
     exit( 1 );
  }
  free( vx );
  free( vc );
  free( vr );
}

void MakeRows( int level, u64 rem, int* vrb, int* vre, 
    u64* vcb, u64* vce, u64* vxb, u64* vxe, int single, u64* countp )
{
  u64 must_bit = ( 1ULL <<  topbit( rem ) );

  for(   ; vrb <  vre - level;  vrb ++ ) {

     u64 row_vec = vecs[ *vrb ];
     if( ! ( row_vec & must_bit ) ) break;

     u64* vcp = vcb;
     u64* vcq = vce;
     while( vcp < vce ) {
        u64 x = row_vec & ( *vcp );
        *vcq = *vcp;
        vcq += ( ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 ) );
        vcp ++;
     }
     if( vcq - vce < N ) continue;

     u64* vxp = vxb;
     u64* vxq = vxe;
     while( vxp < vxe ) {
        u64 x = row_vec & ( *vxp );
        *vxq = *vxp;
        vxq += ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 );
        vxp ++;
     }
     if ( vxq - vxe < 2 ) continue;

     row_vecs[ level ] = row_vec;
     u64 ve = vec_elem[ *vrb ];
     for( int j = 0; j < N; ++ j ) {
        which_row[ ve & 0x3f ] = level;
        ve >>= 6;
     }

     if( level > 1 ) {

        int* vrp = vrb + 1;
        int* vrq = vre;
        while( vrp < vre ) {
           *vrq = *vrp;
           vrq += ( row_vec & ( vecs[ *vrp ] ) ) ? 0 : 1;
           vrp ++;
        }
        // if( vrq - vre < level ) continue; // Don't use continue here!!!

#ifdef V2
        fprintf( stderr,
                " level: %d %09llx  %llu \n", level, row_vec, *countp );
#endif
        MakeRows( level - 1, rem ^ row_vec, vre, vrq, vce, vcq, vxe, vxq, 
                                    single || ( rev_max == row_vec ), countp );
     } else {
        row_vecs[ 0 ] = rem ^ row_vec;
        if( ( ! ( row_vecs[ 0 ] & 1 ) ) || 
              row_vecs[ N - 1 ] >= Reverse( row_vecs[ 0 ] ) ) 
        {  
          sub_count = 0U;
          MakeCols( N - 1, AllBits, vce, vcq, vxe, vxq, 
                    single || rev_max == row_vec || rev_max == row_vecs[ 0 ] ); 
          *countp += sub_count;
        }
     }
     ve = vec_elem[ *vrb ];
     for( int j = 0; j < N; ++ j ) {
        which_row[ ve & 0x3f ] = 0;
        ve >>= 6;
     }
  }
}

void MakeCols( int level, u64 rem,    // 46 / 46
        u64* vcb, u64* vce,  u64* vxb, u64* vxe, int single )
{
  u64 must_bit = ( 1ULL << topbit( rem ) );

  for(    ; vcb < vce - level; vcb ++ ) {

     u64 col_vec = *vcb;
     if( ! ( col_vec & must_bit ) ) break;

     col_vecs[ level ] = col_vec;

     if( level > 1 ) {

        u64* vcp = vcb + 1;  		//  3 / 46
        u64* vcq = vce;
        while( vcp < vce ) { 
           *vcq = *vcp;
           if ( ! ( col_vec & ( *vcp ) ) ) vcq ++ ;
           vcp ++;
        } 
        if( vcq - vce < level ) continue;

        u64* vxp = vxb;              		// 10 / 46
        u64* vxq = vxe;
        while( vxp < vxe ) {                
           u64 x = col_vec & ( *vxp );
           *vxq = *vxp;
           vxq += ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 );
           vxp ++;
        }
        int ndiags = vxq - vxe;
        if ( ndiags < 2 ) continue;

        if( vcq - vce == level ) {
           MakeColsLast( level - 1, rem ^ col_vec, 
                   vce,      vxe, vxq, single || ( rev_max == col_vec ) );
        }
        else {
           MakeCols( level - 1, rem ^ col_vec, 
                   vce, vcq, vxe, vxq, single || ( rev_max == col_vec ) );
        }
     } else {

        int ndiags = DiagPos( vxb, vxe, col_vec, (unsigned*) vxe );
        if ( ndiags < 2 ) continue;

        u64 col_vecs0 = rem ^ col_vec; // col_vecs[ 0 ] is not set.

        if( ( col_vecs0 & 1 ) &&
            row_vecs[ N - 1 ] < Reverse( col_vecs0 ) ) continue;

        int shift = 
           ( single || rev_max == col_vec || rev_max == col_vecs0 ) ? 0 : 1;

        sub_count += MakeDiag( (unsigned*) vxe, ndiags ) << shift;
    }
  }
  return;
}

void MakeColsLast( int level, u64 rem,	//  8 / 46
            u64* cols, u64* vxb, u64* vxe, int single )
{
   u64 rem_inv = ~ rem;
   for( int i = 0; i < level; i ++ ) {
      u64 col_vec = cols[ i ];
      if( rem_inv & col_vec ){ /* printf( "!" ); */ return; }
      rem_inv ^= col_vec;
      col_vecs[ level - i ] = col_vec; 	// col_vecs[0] is not set.
      single = single || ( rev_max == col_vec );
   }
   if( rem_inv & cols[ level ] ){ return; }
   int shift = ( single || ( rev_max == cols[ level ] ) ) ? 0 : 1;

   for( int c = 1; c < level/* + 1*/; c ++ ) { // Don't use col_vecs[0].
      u64* vxq = vxb;
      u64* vxp = vxb;
      while( vxp < vxe ) {
         u64 x = ( *vxp ) & col_vecs[ c ];
         *vxq = *vxp;
         vxq += ( ( x & (x-1) ) ? 0 : 1 ) - ( x ? 0 : 1 );
         vxp ++;
      }
      vxe = vxq;
      if( vxe - vxb < 2 ) return;
   }

   int ndiags = DiagPos( vxb, vxe, col_vecs[ level ], (unsigned*) vxb );
   if( ndiags < 2 ) return;

   sub_count += MakeDiag( (unsigned*) vxb, ndiags ) << shift;
}

int DiagPos( u64* vxb, u64* vxe, u64 col, unsigned* cxb ) 
{
  u64* vxp = vxb;
  unsigned* cxp = cxb;
  while( vxp < vxe ) {                
     u64 x = col & ( *vxp );
     if( ( x & (x-1) ) == 0 && x != 0 ) {
        unsigned pos_code = 0U;
        int posv = 076543210; // octal number
        int onev = 011111111; // octal number
        for( int c = 1; c < N; ++ c ) { // Don't use col_vecs[ 0 ].
           int r = which_row[ ( *vxp & col_vecs[ c ] ) % 37 ];
           int shift = r * OW;
           r = ( ( posv >> shift ) & OMASK );
           shift += OW;
           posv -= ( ( onev >> shift ) << shift );
           pos_code = ( pos_code + r ) * ( N - c );
        }
        *(cxp ++) = pos_code;
     }
     vxp ++;
  }
  return cxp - cxb;
}

unsigned MakeDiag( unsigned* diag_pos, int ndiags )
{
   unsigned  diag_count = 0;
   for( int dx = 0; dx < ndiags - 1; ++ dx ) { 
     unsigned* cx_chk = diag_chk_a[ diag_pos[ dx ] ];
     for( int cx = dx + 1; cx < ndiags; ++ cx ) {
        unsigned cx_pos = diag_pos[ cx ];
        diag_count += ( cx_chk[ cx_pos >> 5 ] >> ( cx_pos & 0x1f ) ) & 1;
     }
   }
   return diag_count;
}

int main( int argc, char* argv[] )
{
   fprintf( stderr, "Size of long long unsigned : %lu \n", sizeof( u64 ) );
   if( sizeof( u64 ) * 8 < NN ) {
      fprintf( stderr, 
               "Size of L L U is too small : %lu \n", sizeof( u64 ) );
      exit( 1 );
   }
   if( check_topbit( ) ) {
      fprintf( stderr, "topbit() does't work correctly on this platform.\n" );
      exit( 1 );
   }

   GenVecs( NN, ( ( NN + 1 ) * N ) / 2, N - 1, 0ULL );
   fprintf( stderr, " nvecs %d\n", nvecs );

   MakeDiagCheck( N );

   u64 row1P = strtoll( argv[ 2 ], NULL, 16 );
   u64 row2P = strtoll( argv[ 4 ], NULL, 16 );
   fprintf( stderr, "r1:%09llx r2:%09llx\n", row1P, row2P );

   u64 total_count = 0ULL;

   for( max_row = 0; max_row < nvecs; ++ max_row ) {
     u64 max_rowP = vecs[ max_row ];
     if( row1P != 0ULL && row1P < max_rowP ) continue;
     if( max_rowP < row1P ) break;
     if( ! ( max_rowP & 1ULL << ( NN - 1 ) ) ) break;
     if( max_rowP < Reverse( vecs[ max_row ] ) ) continue;
     for( second_row = max_row + 1; second_row < nvecs; ++ second_row ) {
         u64 second_rowP = vecs[ second_row ];
         if( row2P != 0ULL && row2P < second_rowP ) continue;
         if( second_rowP < row2P ) break;
         if( max_rowP & second_rowP ) continue;
         if( max_rowP < Reverse( vecs[ second_row ] ) ) continue;
         u64 rem = AllBits ^ max_rowP ^ second_rowP;
         if( rem > second_rowP ) break;

         count_g = 0ULL;
         Count( );
         total_count += count_g;

         printf( "%05d %09llx ", max_row, vecs[ max_row ] );
         printf( "%05d %09llx ", second_row, vecs[ second_row ] );
         printf( "%15llu ", count_g ); 

         #ifndef noMD5

            MD5_CTX md5ctx;
            unsigned char md5bin[ MD5_DIGEST_LENGTH ];
            char md5str[ MD5_DIGEST_LENGTH * 2 + 1 ];

            const char* sample = "e00000007 1c0000038    494199069345";
            int msg_len = strlen( sample );
            char msg[ msg_len + 1 ];
            sprintf( msg, "%09llx %09llx %15llu",
                       vecs[ max_row ], vecs[ second_row ], count_g );

            MD5_Init( &md5ctx );
            MD5_Update( &md5ctx, msg, msg_len );
            MD5_Final( md5bin, &md5ctx );

            for( int i = 0; i < MD5_DIGEST_LENGTH; i ++ )
              sprintf( md5str + ( i * 2 ), "%02x", md5bin[ i ] );
            md5str[ MD5_DIGEST_LENGTH * 2 ] = 0;

            printf( "%s", md5str );

         #endif
         printf( "\n" );
     }
   }
   if( row2P == 0ULL ) fprintf( stderr, "Total: %17llu\n", total_count );

   return 0;
}

</file>