Differences

This shows you the differences between two versions of the page.

--- ms-230905.cu [2023/09/06 22:36] – mino
+++ ms-230905.cu [2023/09/06 22:37] (current) – removed mino
@@ Line 1: / Line 1: @@
-<file c ms-20230131-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;
-}