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-20230131-m.c>
/*

  See https://magicsquare6.net/ .
  
  Updated on 2023/01/31.
  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;
}

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 ======================= */

int diag_chk_size;
u64** diag_chk;

#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( u64* diags, int ndiags );
unsigned MakeDiag2( u64 dx, u64 cx );

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 );

   diag_chk_size = ( perm[ 0 ] >> mask_w ) + 1;
   int v_length = ( perm[ 0 ] >> ( mask_w + 6 ) ) + 1;
   unsigned lower_mask = ( 1 << mask_w ) - 1;
   unsigned pos_mask = ( 1 << ( mask_w * (n-1) ) ) - 1;

   fprintf( stderr, " order %d  diag_chk_size %d  v_length %d \n", 
                                      n, diag_chk_size, v_length );

   diag_chk = (u64**) malloc( sizeof(u64*) * diag_chk_size );
   for( int i = 0; i < diag_chk_size; i ++ ) diag_chk[ i ] = NULL;

   for( int dx = 0; dx < perm_size; dx ++ ) {
      u64* v = (u64*) malloc( sizeof(u64) * v_length );
      for( int j = 0; j < v_length; j ++ ) {
         v[ j ] = 0ULL;
      }
      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 ];
         for( int k = 0; k < n; k ++ ) {
            r[ k ] = ( perm[ cx ] >> ( (n-1-k)*mask_w ) ) & lower_mask;
         }
         if( r[ c[ r[ c[ 0 ] ] ] ] == 0 &&
             r[ c[ r[ c[ 1 ] ] ] ] == 1 &&
             r[ c[ r[ c[ 2 ] ] ] ] == 2 )
         {
            u64 masked = perm[ cx ] & pos_mask;
            v[ masked >> 6 ] |= 1ULL << ( masked & ( ( 1 << 6 ) - 1 ) );
         }
      }
      diag_chk [ perm[ dx ] & pos_mask ] = v; 
   }
}

void FreeDiagCheck( )
{
   for( int i = 0; i < diag_chk_size; i ++ ) {
      if( diag_chk[ i ] ) free( diag_chk[ i ] );
   }
   free( diag_chk );
}

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;

     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( 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;

        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 {

        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 factor = 
           ( single || rev_max == col_vec || rev_max == col_vecs0 ) ? 1 : 2;

        if( ndiags == 2 )  				// 2 / 46
           sub_count += MakeDiag2( *vxe, *( vxe + 1 ) ) * factor;
        else  
           sub_count += MakeDiag( vxe, ndiags ) * factor;
    }
  }
  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; }
   single = single || ( rev_max == cols[ level ] );

   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;
   }
   int ndiags = vxe - vxb;
   if( ndiags < 2 ) return;

   if( ndiags == 2 )  
      sub_count += MakeDiag2( *vxb, *( vxb + 1 ) ) * ( single ? 1 : 2 );
   else  
      sub_count += MakeDiag( vxb, ndiags ) * ( single ? 1 : 2 );
}

unsigned MakeDiag( u64* diags, int ndiags )
{
   unsigned  diag_count = 0;

   unsigned diag_pos[ ndiags ];
   for( int j = 0; j < ndiags; ++ j ) {  // 2 / 22
      unsigned pos_code = 0U;
      for( int c = 1; c < N; ++ c ) { // Don't use col_vecs[ 0 ].
         int r = which_row[ ( diags[ j ] & col_vecs[ c ] ) % 37 ];
         pos_code <<= 3;
         pos_code |= r;
      }
      diag_pos[ j ] = pos_code;
   }

   for( int dx = 0; dx < ndiags - 1; ++ dx ) { 
     u64* cx_chk = diag_chk[ diag_pos[ dx ] ];
     for( int cx = dx + 1; cx < ndiags; ++ cx ) {
        u64 dx_cx = diags[ dx ] & diags[ cx ];
#if N%2
        if( dx_cx && ( ! ( dx_cx & ( dx_cx - 1 ) ) ) ) 
#else
        if( ! dx_cx )
#endif
        {
           unsigned cx_pos = diag_pos[ cx ];
           // printf( " cx %d diag_pos[ cx ]  %d \n", cx, diag_pos[ cx ] );
           diag_count += ( cx_chk[ cx_pos >> 6 ] >> ( cx_pos & 0x3f ) ) & 1;
        }
     }
   }
   return diag_count;
}

unsigned MakeDiag2( u64 dx, u64 cx )
{

   u64 x = dx & cx;
#if N%2
   if( ( ! x ) || ( x & ( x - 1 ) ) ) 
#else
   if( x )
#endif
   return 0;
    
   unsigned dx_pos = 0U;
   unsigned cx_pos = 0U;
   for( int col = 1; col < N; ++ col ) { // Don't use col_vecs[ 0 ].
      int dr = which_row[ ( dx & col_vecs[ col ] ) % 37 ];
      dx_pos <<= 3;
      dx_pos |= dr;
      int cr = which_row[ ( cx & col_vecs[ col ] ) % 37 ];
      cx_pos <<= 3;
      cx_pos |= cr;
   }
   return ( diag_chk[ dx_pos ][ cx_pos >> 6 ] >> ( cx_pos & 0x3f ) ) & 1;
}

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 );

   FreeDiagCheck();
   return 0;
}
</file>