/* to get the segfault compile with
gcc ex.c -lgsl -lgscblas -O2
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <signal.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#define N 16
#define NTOT 76
#define VAL 30
#define SC_EIG 2
#define EPS 0.000001
#define G6LEN(n) (((n)*((n)-1)/2+5)/6+1)
#define BIAS6 63
#define TOPBIT6 32
static gsl_matrix * adj;
static gsl_matrix * partitioned_adj;
static gsl_matrix * adj_ext;
static gsl_vector_complex * eval_ext;
static gsl_matrix_complex * evec_ext;
static gsl_eigen_nonsymmv_workspace * w_ext;
static char gcode[ G6LEN( N) + 3 ] ;
const int eigenvalues[ NTOT] = { 30 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 ,
2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 ,
2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 , - 8 ,
- 8 , - 8 , - 8 , - 8 } ;
char * adjtog6( gsl_matrix * adj, unsigned n) {
unsigned i, j;
int k;
char * p, x;
p = gcode;
* p++ = BIAS6+ n;
k = 6 ;
x = 0 ;
for ( j = 1 ; j < n; ++ j) {
for ( i = 0 ; i < j; ++ i) {
x <<= 1 ;
if ( gsl_matrix_get( adj, j, i) == 1 ) x |= 1 ;
if ( -- k == 0 ) {
* p++ = BIAS6 + x;
k = 6 ;
x = 0 ;
}
}
}
if ( k != 6 ) * p++ = BIAS6 + ( x << k) ;
* p++ = '\n ' ;
* p = '\0 ' ;
return gcode;
}
static gsl_matrix * partitioned_am( gsl_matrix * adj, gsl_matrix * partitioned_adj, unsigned n) {
unsigned i, j;
unsigned total_edges = 0 ;
gsl_matrix_set_zero( partitioned_adj) ;
for ( i= 0 ; i < n; i++ ) {
unsigned deg = 0 ;
for ( j = 0 ; j < n; j++ ) {
unsigned el = gsl_matrix_get( adj, i, j) ;
if ( el) {
gsl_matrix_set( partitioned_adj, i, j, 1 ) ;
deg+= 1 ;
}
}
gsl_matrix_set( partitioned_adj, i, n, VAL- deg) ;
gsl_matrix_set( partitioned_adj, n, i, ( double ) ( VAL- deg) / ( NTOT- n) ) ;
total_edges += deg;
}
gsl_matrix_set( partitioned_adj, n, n, ( double ) 2 * ( NTOT* VAL/ 2 + total_edges/ 2 - n* VAL) / ( NTOT- n) ) ;
return partitioned_adj;
}
static int cmp( const void * elem1, const void * elem2) {
const double * a = elem1, * b = elem2;
return * a > * b ? - 1 : * a < * b ? 1 : 0 ;
}
static double * spectrum( gsl_matrix * m, gsl_vector_complex * eval, gsl_matrix_complex * evec,
gsl_eigen_nonsymmv_workspace * w, unsigned size) {
static double eigs[ N+ 2 ] ;
unsigned i;
gsl_eigen_nonsymmv ( m, eval, evec, w) ;
for ( i = 0 ; i < size ; i++ ) {
gsl_complex eval_i = gsl_vector_complex_get ( eval, i) ;
eigs[ i] = GSL_REAL( eval_i) ;
}
qsort ( eigs
, size
, sizeof ( double ) , cmp
) ;
return eigs;
}
static unsigned does_interlace( double eigs_sc[ N+ 2 ] , unsigned n) {
unsigned i;
for ( i = 0 ; i < n; i++ ) {
double expr = eigenvalues[ i] - eigs_sc[ i] ;
if ( expr
<= EPS
&& fabs ( expr
) >= EPS
) { return 0 ;
}
expr = eigs_sc[ i] - eigenvalues[ NTOT- n+ i] ;
if ( expr
<= EPS
&& fabs ( expr
) >= EPS
) { return 0 ;
}
}
return 1 ;
}
static unsigned isInterlaced( gsl_matrix * adj) {
gsl_matrix * m = partitioned_am( adj, partitioned_adj, N+ 1 ) ;
double * eigs = spectrum( m, eval_ext, evec_ext, w_ext, N+ 2 ) ;
return does_interlace( eigs, N+ 2 ) ;
}
static void expand( gsl_matrix * adj, unsigned n) {
unsigned i, j;
gsl_matrix * adj_ext = gsl_matrix_alloc( n+ 1 , n+ 1 ) ;
gsl_matrix_set_zero( adj_ext) ;
for ( i = 0 ; i < n; i++ ) {
for ( j = i+ 1 ; j < n; j++ ) {
double val = gsl_matrix_get( adj, i, j) ;
if ( val) {
gsl_matrix_set( adj_ext, i, j, 1.0 ) ;
gsl_matrix_set( adj_ext, j, i, 1.0 ) ;
}
}
}
for ( i = 0 ; i < ( 1U << n) ; i++ ) {
for ( j = 0 ; j < n; j++ ) {
if ( i & ( 1U << j) ) {
gsl_matrix_set( adj_ext, n, j, 1 ) ;
gsl_matrix_set( adj_ext, j, n, 1 ) ;
}
}
if ( isInterlaced( adj_ext) ) {
printf ( "%s" , adjtog6
( adj_ext
, n
+ 1 ) ) ; }
for ( j = 0 ; j < n; j++ ) {
if ( i & ( 1 << j) ) {
gsl_matrix_set( adj_ext, n, j, 0 ) ;
gsl_matrix_set( adj_ext, j, n, 0 ) ;
}
}
}
gsl_matrix_free( adj_ext) ;
}
static void stringtoadj( char * s) {
char * p;
int i, j, k, x = 0 ;
gsl_matrix_set_zero( adj) ;
p = s + 1 ;
k = 1 ;
for ( j = 1 ; j < N; ++ j) {
for ( i = 0 ; i < j; ++ i) {
if ( -- k == 0 ) {
k = 6 ;
x = * ( p++ ) - BIAS6;
}
if ( x & TOPBIT6) {
gsl_matrix_set( adj, i, j, 1 ) ;
gsl_matrix_set( adj, j, i, 1 ) ;
}
x <<= 1 ;
}
}
}
int main( int argc, char * argv[ ] ) {
adj = gsl_matrix_calloc( N, N) ;
partitioned_adj = gsl_matrix_calloc( N+ 2 , N+ 2 ) ;
adj_ext = gsl_matrix_alloc( N+ 1 , N+ 1 ) ;
eval_ext = gsl_vector_complex_alloc( N+ 2 ) ;
evec_ext = gsl_matrix_complex_alloc ( N+ 2 , N+ 2 ) ;
w_ext = gsl_eigen_nonsymmv_alloc( N+ 2 ) ;
if ( ! w_ext || ! evec_ext || ! eval_ext ||! adj_ext) {
return 1 ;
}
if ( ! adj || ! partitioned_adj) {
puts ( "Some error somewhere" ) ; return 1 ;
}
stringtoadj( "O??O@BI@?Ej@hBAt}YFwa" ) ;
expand( adj, N) ;
return 0 ;
}
LyogdG8gZ2V0IHRoZSBzZWdmYXVsdCBjb21waWxlIHdpdGggCgogICBnY2MgZXguYyAtbGdzbCAtbGdzY2JsYXMgLU8yCgoqLyAgIAoKI2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPHN0cmluZy5oPgojaW5jbHVkZSA8c2lnbmFsLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojaW5jbHVkZSA8Z3NsL2dzbF9saW5hbGcuaD4KI2luY2x1ZGUgPGdzbC9nc2xfbWF0cml4Lmg+CiNpbmNsdWRlIDxnc2wvZ3NsX2VpZ2VuLmg+CgojZGVmaW5lIE4gMTYKCiNkZWZpbmUgTlRPVCA3NiAKI2RlZmluZSBWQUwgMzAKCiNkZWZpbmUgU0NfRUlHIDIKCiNkZWZpbmUgRVBTIDAuMDAwMDAxCgojZGVmaW5lIEc2TEVOKG4pICAoKChuKSooKG4pLTEpLzIrNSkvNisxKSAKI2RlZmluZSBCSUFTNiA2MwojZGVmaW5lIFRPUEJJVDYgMzIKCnN0YXRpYyBnc2xfbWF0cml4ICphZGo7CnN0YXRpYyBnc2xfbWF0cml4ICpwYXJ0aXRpb25lZF9hZGo7CgpzdGF0aWMgZ3NsX21hdHJpeCAqYWRqX2V4dDsKCnN0YXRpYyBnc2xfdmVjdG9yX2NvbXBsZXggKmV2YWxfZXh0OwpzdGF0aWMgZ3NsX21hdHJpeF9jb21wbGV4ICpldmVjX2V4dDsKc3RhdGljIGdzbF9laWdlbl9ub25zeW1tdl93b3Jrc3BhY2UgKndfZXh0OwoKc3RhdGljIGNoYXIgZ2NvZGVbRzZMRU4oTikrM107IAoKY29uc3QgaW50IGVpZ2VudmFsdWVzW05UT1RdID0gezMwLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAyLCAKICAgIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIDIsIAogICAgMiwgMiwgMiwgMiwgMiwgMiwgMiwgMiwgMiwgMiwgLTgsIC04LCAtOCwgLTgsIC04LCAtOCwgLTgsIC04LCAtOCwgLTgsIC04LCAtOCwgLTgsIC04LCAKICAgIC04LCAtOCwgLTgsIC04fTsKCgpjaGFyICphZGp0b2c2KGdzbF9tYXRyaXggKmFkaiwgdW5zaWduZWQgbikgewoKICAgIHVuc2lnbmVkIGksajsKICAgIGludCAgazsKICAgIAogICAgY2hhciAqcCx4OwoKICAgIHAgPSBnY29kZTsKICAgICpwKysgPSBCSUFTNituOwoKICAgIGsgPSA2OwogICAgeCA9IDA7CiAgICAKICAgIGZvciAoaiA9IDE7IGogPCBuOyArK2opIHsKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgajsgKytpKSB7CiAgICAgICAgICAgIHggPDw9IDE7CiAgICAgICAgICAgIGlmIChnc2xfbWF0cml4X2dldChhZGosIGosaSkgPT0gMSkgeCB8PSAxOwoKICAgICAgICAgICAgaWYgKC0tayA9PSAwKSB7CiAgICAgICAgICAgICAgICAqcCsrID0gQklBUzYgKyB4OwogICAgICAgICAgICAgICAgayA9IDY7CiAgICAgICAgICAgICAgICB4ID0gMDsKICAgICAgICAgICAgfQogICAgICAgIH0KICAgIH0KCiAgICBpZiAoayAhPSA2KSAqcCsrID0gQklBUzYgKyAoeCA8PCBrKTsKCiAgICAqcCsrID0gJ1xuJzsKICAgICpwID0gJ1wwJzsKCiAgICByZXR1cm4gZ2NvZGU7Cn0KCnN0YXRpYyBnc2xfbWF0cml4ICpwYXJ0aXRpb25lZF9hbShnc2xfbWF0cml4ICphZGosIGdzbF9tYXRyaXggKnBhcnRpdGlvbmVkX2FkaiwgdW5zaWduZWQgbikgewoKICAgIHVuc2lnbmVkIGksajsKCiAgICB1bnNpZ25lZCB0b3RhbF9lZGdlcyA9IDA7IAoKICAgIGdzbF9tYXRyaXhfc2V0X3plcm8ocGFydGl0aW9uZWRfYWRqKTsKCiAgICBmb3IgKGk9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICB1bnNpZ25lZCBkZWcgPSAwOwogICAgICAgIGZvciAoaiA9IDA7IGogPCBuOyBqKyspIHsKICAgICAgICAgICAgdW5zaWduZWQgZWwgPSBnc2xfbWF0cml4X2dldChhZGosaSxqKTsKICAgICAgICAgICAgaWYgKGVsKSAgewogICAgICAgICAgICAgICAgZ3NsX21hdHJpeF9zZXQocGFydGl0aW9uZWRfYWRqLCBpLGosMSk7CiAgICAgICAgICAgICAgICBkZWcrPTE7CiAgICAgICAgICAgIH0gICAKICAgICAgICB9CiAgICAgICAgZ3NsX21hdHJpeF9zZXQocGFydGl0aW9uZWRfYWRqLGksIG4sIFZBTC1kZWcpOwogICAgICAgIGdzbF9tYXRyaXhfc2V0KHBhcnRpdGlvbmVkX2FkaiwgbixpLCAoZG91YmxlKShWQUwtZGVnKS8oTlRPVC1uKSk7CiAgICAgICAgdG90YWxfZWRnZXMgKz0gZGVnOwogICAgfQogICAgZ3NsX21hdHJpeF9zZXQocGFydGl0aW9uZWRfYWRqLG4sIG4sIChkb3VibGUpIDIqKE5UT1QqVkFMLzIgKyB0b3RhbF9lZGdlcy8yIC0gbipWQUwpLyhOVE9ULW4pKTsKCiAgICByZXR1cm4gcGFydGl0aW9uZWRfYWRqOwp9CgpzdGF0aWMgaW50IGNtcChjb25zdCB2b2lkKiBlbGVtMSwgY29uc3Qgdm9pZCogZWxlbTIpIHsKCiAgICBjb25zdCBkb3VibGUgKmEgPSBlbGVtMSwgKmIgPSBlbGVtMjsgCiAgICAKICAgIHJldHVybiAqYSA+ICpiID8gLTEgOiAqYSA8ICpiID8gMSA6IDA7Cn0KCgpzdGF0aWMgZG91YmxlICpzcGVjdHJ1bShnc2xfbWF0cml4ICptLCBnc2xfdmVjdG9yX2NvbXBsZXggKmV2YWwsIGdzbF9tYXRyaXhfY29tcGxleCAqZXZlYywgCiAgICAgICAgZ3NsX2VpZ2VuX25vbnN5bW12X3dvcmtzcGFjZSAqdywgdW5zaWduZWQgc2l6ZSkgewoKICAgIHN0YXRpYyBkb3VibGUgZWlnc1tOKzJdOwogICAgdW5zaWduZWQgaTsKCiAgICBnc2xfZWlnZW5fbm9uc3ltbXYgKG0sIGV2YWwsIGV2ZWMsIHcpOyAgICAKCiAgICBmb3IoaSA9IDA7IGkgPCBzaXplIDsgaSsrKSB7CiAgICAgICAgZ3NsX2NvbXBsZXggZXZhbF9pID0gZ3NsX3ZlY3Rvcl9jb21wbGV4X2dldCAoZXZhbCwgaSk7CiAgICAgICAgZWlnc1tpXSA9IEdTTF9SRUFMKGV2YWxfaSk7CiAgICB9CgogICAgcXNvcnQoZWlncyxzaXplLCBzaXplb2YoZG91YmxlKSxjbXApOyAKICAgIAogICAgcmV0dXJuIGVpZ3M7Cn0KCnN0YXRpYyB1bnNpZ25lZCBkb2VzX2ludGVybGFjZShkb3VibGUgZWlnc19zY1tOKzJdLCB1bnNpZ25lZCBuKSB7CgogICAgdW5zaWduZWQgaTsKCiAgICBmb3IgKGkgPSAwOyBpIDwgbjsgaSsrKSB7CgogICAgICAgIGRvdWJsZSBleHByID0gZWlnZW52YWx1ZXNbaV0tZWlnc19zY1tpXTsKCiAgICAgICAgaWYgKGV4cHIgPD0gRVBTICYmIGZhYnMoZXhwcikgPj0gRVBTKSB7IAogICAgICAgICAgICByZXR1cm4gMDsgICAgCiAgICAgICAgfQoKICAgICAgICBleHByID0gZWlnc19zY1tpXSAtIGVpZ2VudmFsdWVzW05UT1QtbitpXTsKCiAgICAgICAgaWYgKGV4cHIgPD0gRVBTICYmIGZhYnMoZXhwcikgPj0gRVBTKSB7CiAgICAgICAgICAgIHJldHVybiAwOwogICAgICAgIH0gCiAgICB9CgogICAgcmV0dXJuIDE7IAp9CgoKc3RhdGljIHVuc2lnbmVkIGlzSW50ZXJsYWNlZChnc2xfbWF0cml4ICphZGopIHsKCiAgICBnc2xfbWF0cml4ICptID0gcGFydGl0aW9uZWRfYW0oYWRqLHBhcnRpdGlvbmVkX2FkaixOKzEpOwoKICAgIGRvdWJsZSAqZWlncyA9IHNwZWN0cnVtKG0sIGV2YWxfZXh0LCBldmVjX2V4dCwgd19leHQsIE4rMik7CiAgICAKICAgIHJldHVybiBkb2VzX2ludGVybGFjZShlaWdzLE4rMik7Cn0KCnN0YXRpYyB2b2lkIGV4cGFuZChnc2xfbWF0cml4ICphZGosIHVuc2lnbmVkIG4pIHsKCiAgICB1bnNpZ25lZCBpLGo7CgogICAgZ3NsX21hdHJpeCAqYWRqX2V4dCA9IGdzbF9tYXRyaXhfYWxsb2MobisxLG4rMSk7CiAgICBnc2xfbWF0cml4X3NldF96ZXJvKGFkal9leHQpOwoKICAgIGZvciAoaSA9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICBmb3IgKGogPSBpKzE7IGogPCBuOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIHZhbCA9IGdzbF9tYXRyaXhfZ2V0KGFkaixpLGopOwogICAgICAgICAgICBpZiAodmFsKSB7CiAgICAgICAgICAgICAgICBnc2xfbWF0cml4X3NldChhZGpfZXh0LCBpLGosIDEuMCk7CiAgICAgICAgICAgICAgICBnc2xfbWF0cml4X3NldChhZGpfZXh0LCBqLGksIDEuMCk7CiAgICAgICAgICAgIH0KICAgICAgICB9CiAgICB9CiAgICAKICAgIGZvciAoaSA9IDA7IGkgPCAoMVUgPDxuKTsgaSsrKSB7CgogICAgICAgIGZvciAoaiA9IDA7IGogPCBuOyBqKyspIHsKICAgICAgICAgICAgaWYgKGkgJiAoMVUgPDwgaikgKSB7CiAgICAgICAgICAgICAgICBnc2xfbWF0cml4X3NldChhZGpfZXh0LCBuLCBqLDEpOwogICAgICAgICAgICAgICAgZ3NsX21hdHJpeF9zZXQoYWRqX2V4dCwgaiwgbiwxKTsKICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAgICAgaWYgKGlzSW50ZXJsYWNlZChhZGpfZXh0KSkgewogICAgICAgICAgICBwcmludGYoIiVzIiwgYWRqdG9nNihhZGpfZXh0LCBuKzEpKTsKICAgICAgICB9CgogICAgICAgIGZvciAoaiA9IDA7IGogPCBuOyBqKyspIHsKICAgICAgICAgICAgaWYgKGkgJiAoMSA8PCBqKSApIHsKICAgICAgICAgICAgICAgIGdzbF9tYXRyaXhfc2V0KGFkal9leHQsIG4sIGosMCk7CiAgICAgICAgICAgICAgICBnc2xfbWF0cml4X3NldChhZGpfZXh0LCBqLCBuLDApOwogICAgICAgICAgICB9CiAgICAgICAgfQoKICAgIH0KICAgIGdzbF9tYXRyaXhfZnJlZShhZGpfZXh0KTsKfQoKc3RhdGljIHZvaWQgc3RyaW5ndG9hZGooY2hhciAqcykgewoKCWNoYXIgKnA7CglpbnQgaSxqLGsseCA9IDA7CgogICAgZ3NsX21hdHJpeF9zZXRfemVybyhhZGopOwoKICAgIHAgPSBzICsgMTsKCiAgICBrID0gMTsKICAgIGZvciAoaiA9IDE7IGogPCBOOyArK2opIHsKICAgICAgICBmb3IgKGkgPSAwOyBpIDwgajsgKytpKSB7CiAgICAgICAgICAgIGlmICgtLWsgPT0gMCkgewoJCSAgICAgICAgayA9IDY7CgkJICAgICAgICB4ID0gKihwKyspIC0gQklBUzY7CiAgICAgICAgICAgIH0KCSAgICAKICAgICAgICAgICAgaWYgKHggJiBUT1BCSVQ2KSB7CiAgICAgICAgICAgICAgICBnc2xfbWF0cml4X3NldChhZGosaSxqLDEpOwogICAgICAgICAgICAgICAgZ3NsX21hdHJpeF9zZXQoYWRqLGosaSwxKTsKICAgICAgICAgICAgfQogICAgICAgICAgICB4IDw8PSAxOwogICAgICAgIH0KICAgIH0KfQoKaW50IG1haW4oaW50IGFyZ2MsIGNoYXIgKmFyZ3ZbXSkgewoKICAgIGFkaiA9IGdzbF9tYXRyaXhfY2FsbG9jKE4sIE4pOwogICAgcGFydGl0aW9uZWRfYWRqID0gZ3NsX21hdHJpeF9jYWxsb2MoTisyLCBOKzIpOwoKICAgIGFkal9leHQgPSBnc2xfbWF0cml4X2FsbG9jKE4rMSwgTisxKTsKCiAgICBldmFsX2V4dCA9IGdzbF92ZWN0b3JfY29tcGxleF9hbGxvYyhOKzIpOwogICAgZXZlY19leHQgPSBnc2xfbWF0cml4X2NvbXBsZXhfYWxsb2MgKE4rMiwgTisyKTsKICAgIHdfZXh0ID0gZ3NsX2VpZ2VuX25vbnN5bW12X2FsbG9jKE4rMik7CiAgICAKICAgIGlmICghd19leHQgfHwgIWV2ZWNfZXh0IHx8ICFldmFsX2V4dCB8fCFhZGpfZXh0KSB7CiAgICAgICAgcHV0cygiRXJycnIiKTsKICAgICAgICByZXR1cm4gMTsKICAgIH0gICAgCgogICAgaWYgKCFhZGogfHwgIXBhcnRpdGlvbmVkX2FkaikgeyAKICAgICAgICBwdXRzKCJTb21lIGVycm9yIHNvbWV3aGVyZSIpOwogICAgICAgIHJldHVybiAxOwogICAgfQogICAgCiAgICBzdHJpbmd0b2FkaigiTz8/T0BCSUA/RWpAaEJBdH1ZRndhIik7CiAgICBleHBhbmQoYWRqLE4pOwoKICAgIHJldHVybiAwOwoKfQ==