2016-11-20 22:19:08 +00:00
# include "weyl.h"
# include "queue.h"
# include <stdio.h>
# include <memory.h>
# include <stdlib.h>
# define BIT(n) ((uint64_t)1 << (n))
typedef struct {
weylid_t id ;
int position ;
} weylid_lookup_t ;
static int search ( const void * key , const void * base , size_t nmem , size_t size , int ( * compar ) ( const void * , const void * , void * ) , void * arg ) ;
static int compare_root_vectors ( int rank , const int * x , const int * y ) ;
static int compare_root_vectors_qsort ( const void * x , const void * y , void * arg ) ;
static int compare_weylid_lookup ( const void * x , const void * y ) ;
static int lookup_id ( weylid_t id , weylid_lookup_t * list , int len ) ;
static weylid_t multiply_generator ( int s , weylid_t w , const int * simple , const int * mapping , int rank , int positive ) ;
static void reflect_root_vector ( const int * cartan , int rank , int i , int * old , int * new ) ;
/***************** simple helper functions **********************************/
// glibc search function, but with user pointer and returning index (or -1 if not found)
static int search ( const void * key , const void * base , size_t nmemb , size_t size , int ( * compar ) ( const void * , const void * , void * ) , void * arg )
{
size_t l , u , idx ;
const void * p ;
int comparison ;
l = 0 ;
u = nmemb ;
while ( l < u ) {
idx = ( l + u ) / 2 ;
p = ( void * ) ( ( ( const char * ) base ) + ( idx * size ) ) ;
comparison = ( * compar ) ( key , p , arg ) ;
if ( comparison < 0 )
u = idx ;
else if ( comparison > 0 )
l = idx + 1 ;
else
return idx ;
}
return - 1 ;
}
// maybe we want a different ordering here?
static int compare_root_vectors ( int rank , const int * x , const int * y )
{
for ( int i = 0 ; i < rank ; i + + )
if ( x [ i ] ! = y [ i ] )
return x [ i ] - y [ i ] ;
return 0 ;
}
static int compare_root_vectors_qsort ( const void * x , const void * y , void * arg )
{
return compare_root_vectors ( * ( ( int * ) arg ) , x , y ) ;
}
static int compare_weylid ( const void * x , const void * y )
{
weylid_t u = * ( ( weylid_t * ) x ) ;
weylid_t v = * ( ( weylid_t * ) y ) ;
return u > v ? 1 : u < v ? - 1 : 0 ;
}
static int compare_weylid_lookup ( const void * x , const void * y )
{
weylid_t u = ( ( weylid_lookup_t * ) x ) - > id ;
weylid_t v = ( ( weylid_lookup_t * ) y ) - > id ;
return u > v ? 1 : u < v ? - 1 : 0 ;
}
static int lookup_id ( weylid_t id , weylid_lookup_t * list , int len )
{
weylid_lookup_t key ;
key . id = id ;
weylid_lookup_t * p = ( weylid_lookup_t * ) bsearch ( & key , list , len , sizeof ( weylid_lookup_t ) , compare_weylid_lookup ) ;
return p - > position ;
}
static weylid_t multiply_generator ( int s , weylid_t w , const int * simple , const int * mapping , int rank , int positive )
{
weylid_t sw = 0 ;
for ( int i = 0 ; i < positive ; i + + ) {
if ( w & BIT ( i ) )
if ( mapping [ i * rank + s ] ! = - 1 )
sw | = BIT ( mapping [ i * rank + s ] ) ;
}
if ( w & BIT ( simple [ s ] ) )
return sw ;
else
return sw | BIT ( simple [ s ] ) ;
}
static void reflect_root_vector ( const int * cartan , int rank , int i , int * old , int * new )
{
memcpy ( new , old , rank * sizeof ( int ) ) ;
for ( int j = 0 ; j < rank ; j + + )
new [ i ] - = cartan [ i * rank + j ] * old [ j ] ;
}
/************* Weyl group infos ************************/
int weyl_rank ( semisimple_type_t type )
{
int rank = 0 ;
for ( int i = 0 ; i < type . n ; i + + )
rank + = type . factors [ i ] . rank ;
return rank ;
}
int weyl_order ( semisimple_type_t type )
{
int order = 1 ;
for ( int i = 0 ; i < type . n ; i + + ) {
switch ( type . factors [ i ] . series ) {
case ' A ' :
for ( int j = 1 ; j < = type . factors [ i ] . rank + 1 ; j + + )
order * = j ;
break ;
case ' B ' : case ' C ' :
for ( int j = 1 ; j < = type . factors [ i ] . rank ; j + + )
order * = 2 * j ;
break ;
case ' D ' :
for ( int j = 2 ; j < = type . factors [ i ] . rank ; j + + )
order * = 2 * j ;
break ;
case ' E ' :
if ( type . factors [ i ] . rank = = 6 )
order * = 51840 ;
else if ( type . factors [ i ] . rank = = 7 )
order * = 2903040 ;
else if ( type . factors [ i ] . rank = = 8 )
order * = 696729600 ;
else
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
break ;
case ' F ' :
ERROR ( type . factors [ i ] . rank ! = 4 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
order * = 1152 ;
break ;
case ' G ' :
ERROR ( type . factors [ i ] . rank ! = 2 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
order * = 12 ;
break ;
default :
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
}
}
return order ;
}
int weyl_hyperplanes ( semisimple_type_t type )
{
int hyperplanes = 0 ;
for ( int i = 0 ; i < type . n ; i + + ) {
switch ( type . factors [ i ] . series ) {
case ' A ' :
hyperplanes + = ( type . factors [ i ] . rank * ( type . factors [ i ] . rank + 1 ) ) / 2 ;
break ;
case ' B ' : case ' C ' :
hyperplanes + = type . factors [ i ] . rank * type . factors [ i ] . rank ;
break ;
case ' D ' :
hyperplanes + = type . factors [ i ] . rank * ( type . factors [ i ] . rank - 1 ) ;
break ;
case ' E ' :
if ( type . factors [ i ] . rank = = 6 )
hyperplanes + = 36 ;
else if ( type . factors [ i ] . rank = = 7 )
hyperplanes + = 63 ;
else if ( type . factors [ i ] . rank = = 8 )
hyperplanes + = 120 ;
else
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
break ;
case ' F ' :
ERROR ( type . factors [ i ] . rank ! = 4 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
hyperplanes + = 24 ;
break ;
case ' G ' :
ERROR ( type . factors [ i ] . rank ! = 2 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
hyperplanes + = 6 ;
break ;
default :
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ i ] . series , type . factors [ i ] . rank ) ;
}
}
return hyperplanes ;
}
int weyl_opposition ( semisimple_type_t type , int simple_root )
{
int offset = 0 ;
int factor = 0 ;
int r , iota_r ;
for ( factor = 0 ; factor < type . n ; factor + + )
if ( simple_root < offset + type . factors [ factor ] . rank )
break ;
else
offset + = type . factors [ factor ] . rank ;
r = simple_root - offset ;
switch ( type . factors [ factor ] . series ) {
2016-12-13 20:26:48 +00:00
case ' B ' : case ' C ' : case ' F ' : case ' G ' :
iota_r = r ;
2016-11-20 22:19:08 +00:00
break ;
2016-12-13 20:26:48 +00:00
case ' A ' :
iota_r = type . factors [ factor ] . rank - 1 - r ;
2016-11-20 22:19:08 +00:00
break ;
case ' D ' :
2016-12-13 20:26:48 +00:00
if ( type . factors [ factor ] . rank % 2 = = 0 )
iota_r = r ;
else
iota_r = r = = 0 ? 1 : r = = 1 ? 0 : r ;
2016-11-20 22:19:08 +00:00
break ;
case ' E ' :
2016-12-13 20:26:48 +00:00
if ( type . factors [ factor ] . rank = = 6 )
iota_r = r ;
else
iota_r = r = = 3 | | r = = 4 ? r : 5 - r ;
2016-11-20 22:19:08 +00:00
break ;
default :
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ factor ] . series , type . factors [ factor ] . rank ) ;
}
return iota_r + offset ;
}
void weyl_cartan_matrix ( semisimple_type_t type , int * m )
{
int offset = 0 ;
int rank = weyl_rank ( type ) ;
int * * A = ( int * * ) malloc ( rank * sizeof ( int * ) ) ;
memset ( m , 0 , rank * rank * sizeof ( int ) ) ;
for ( int i = 0 ; i < rank ; i + + )
m [ i * rank + i ] = 2 ;
for ( int k = 0 ; k < type . n ; k + + ) {
for ( int i = 0 ; i < type . factors [ k ] . rank ; i + + ) // A is the submatrix corresponding to the current simple factor
A [ i ] = & m [ ( i + offset ) * rank + offset ] ;
2016-12-13 20:26:48 +00:00
for ( int i = 1 ; i < type . factors [ k ] . rank ; i + + ) {
A [ i ] [ i - 1 ] = - 1 ;
A [ i - 1 ] [ i ] = - 1 ;
}
2016-11-20 22:19:08 +00:00
switch ( type . factors [ k ] . series ) {
case ' A ' :
break ;
2016-12-13 20:26:48 +00:00
case ' B ' :
ERROR ( type . factors [ k ] . rank < 2 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 0 ] [ 1 ] = - 2 ;
2016-11-20 22:19:08 +00:00
break ;
case ' C ' :
2016-12-13 20:26:48 +00:00
ERROR ( type . factors [ k ] . rank < 2 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 1 ] [ 0 ] = - 2 ;
2016-11-20 22:19:08 +00:00
break ;
case ' D ' :
2016-12-13 20:26:48 +00:00
ERROR ( type . factors [ k ] . rank < 3 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 0 ] [ 1 ] = A [ 1 ] [ 0 ] = 0 ;
A [ 0 ] [ 2 ] = A [ 2 ] [ 0 ] = - 1 ;
2016-11-20 22:19:08 +00:00
break ;
case ' E ' :
2016-12-13 20:26:48 +00:00
ERROR ( type . factors [ k ] . rank < 6 | | type . factors [ k ] . rank > 8 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 1 ] [ 2 ] = A [ 2 ] [ 1 ] = 0 ;
A [ 1 ] [ 3 ] = A [ 3 ] [ 1 ] = - 1 ;
2016-11-20 22:19:08 +00:00
break ;
case ' F ' :
2016-12-13 20:26:48 +00:00
ERROR ( type . factors [ k ] . rank ! = 4 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 3 ] [ 2 ] = - 2 ;
2016-11-20 22:19:08 +00:00
break ;
case ' G ' :
2016-12-13 20:26:48 +00:00
ERROR ( type . factors [ k ] . rank ! = 2 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
A [ 1 ] [ 0 ] = - 3 ;
2016-11-20 22:19:08 +00:00
break ;
default :
ERROR ( 1 , " A Weyl group of type %c%d does not exist or is not implemented! \n " , type . factors [ k ] . series , type . factors [ k ] . rank ) ;
}
offset + = type . factors [ k ] . rank ;
}
free ( A ) ;
}
/************ memory allocation ********************/
weylgroup_element_t * weyl_alloc ( semisimple_type_t type )
{
int rank = weyl_rank ( type ) ;
int order = weyl_order ( type ) ;
int * left = ( int * ) malloc ( rank * order * sizeof ( int ) ) ;
int * right = ( int * ) malloc ( rank * order * sizeof ( int ) ) ;
weylgroup_element_t * group = ( weylgroup_element_t * ) malloc ( order * sizeof ( weylgroup_element_t ) ) ;
for ( int i = 0 ; i < order ; i + + ) {
group [ i ] . left = & left [ i * rank ] ;
group [ i ] . right = & right [ i * rank ] ;
}
return group ;
}
void weyl_free ( weylgroup_element_t * x )
{
free ( x [ 0 ] . left ) ;
free ( x [ 0 ] . right ) ;
free ( x ) ;
}
void weyl_generate ( semisimple_type_t type , weylgroup_element_t * group )
{
int rank , order , positive ;
queue_t queue ;
int current ;
int roots_known , elements , length_elements , nextids_count ;
int * cartan_matrix ;
int * root_vectors ;
int * vector ;
int * simple_roots ;
int * root_mapping ;
weylid_t * ids , * edges , * nextids ;
weylid_lookup_t * lookup ;
rank = weyl_rank ( type ) ;
order = weyl_order ( type ) ;
positive = weyl_hyperplanes ( type ) ;
ERROR ( positive > 64 , " We can't handle root systems with more than 64 positive roots! \n " ) ;
cartan_matrix = ( int * ) malloc ( rank * rank * sizeof ( int ) ) ;
root_vectors = ( int * ) malloc ( 2 * positive * rank * sizeof ( int ) ) ;
vector = ( int * ) malloc ( rank * sizeof ( int ) ) ;
root_mapping = ( int * ) malloc ( positive * rank * sizeof ( int ) ) ;
simple_roots = ( int * ) malloc ( rank * sizeof ( int ) ) ;
ids = ( weylid_t * ) malloc ( order * sizeof ( weylid_t ) ) ;
edges = ( weylid_t * ) malloc ( rank * order * sizeof ( weylid_t ) ) ;
nextids = ( weylid_t * ) malloc ( rank * order * sizeof ( weylid_t ) ) ;
lookup = ( weylid_lookup_t * ) malloc ( order * sizeof ( weylid_lookup_t ) ) ;
weyl_cartan_matrix ( type , cartan_matrix ) ;
// enumerate roots
memset ( root_vectors , 0 , 2 * positive * rank * sizeof ( int ) ) ;
// first the simple roots
queue_init ( & queue ) ;
for ( int i = 0 ; i < rank ; i + + ) {
root_vectors [ rank * i + i ] = 1 ;
queue_put ( & queue , i ) ;
}
// and then we get all others by reflecting
roots_known = rank ;
while ( ( current = queue_get ( & queue ) ) ! = - 1 ) {
for ( int i = 0 ; i < rank ; i + + ) {
reflect_root_vector ( cartan_matrix , rank , i , & root_vectors [ rank * current ] , vector ) ;
int j ;
for ( j = 0 ; j < roots_known ; j + + )
if ( compare_root_vectors ( rank , & root_vectors [ rank * j ] , vector ) = = 0 )
break ;
if ( j = = roots_known ) {
memcpy ( & root_vectors [ rank * roots_known ] , vector , rank * sizeof ( int ) ) ;
queue_put ( & queue , roots_known ) ;
roots_known + + ;
}
}
}
ERROR ( roots_known ! = 2 * positive , " Number of roots does not match! \n " )
// sort roots and restrict to positives
qsort_r ( root_vectors , 2 * positive , rank * sizeof ( int ) , compare_root_vectors_qsort , & rank ) ;
memcpy ( root_vectors , & root_vectors [ positive * rank ] , positive * rank * sizeof ( int ) ) ;
for ( int i = 0 ; i < positive ; i + + ) {
for ( int j = 0 ; j < rank ; j + + ) {
reflect_root_vector ( cartan_matrix , rank , j , & root_vectors [ rank * i ] , vector ) ;
root_mapping [ i * rank + j ] =
search ( vector , root_vectors , positive , rank * sizeof ( int ) , compare_root_vectors_qsort , & rank ) ;
}
}
// where in the list are the simple roots?
for ( int i = 0 ; i < rank ; i + + ) {
memset ( vector , 0 , rank * sizeof ( int ) ) ;
vector [ i ] = 1 ;
simple_roots [ i ] = search ( vector , root_vectors , positive , rank * sizeof ( int ) , compare_root_vectors_qsort , & rank ) ;
}
// enumerate weyl group elements using difference sets
nextids [ 0 ] = 0 ;
nextids_count = 1 ;
elements = 0 ;
for ( int len = 0 ; len < = positive ; len + + ) {
length_elements = 0 ;
// find unique ids in edges added in the last iteration
qsort ( nextids , nextids_count , sizeof ( weylid_t ) , compare_weylid ) ;
for ( int i = 0 ; i < nextids_count ; i + + )
if ( i = = 0 | | nextids [ i ] ! = nextids [ i - 1 ] )
ids [ elements + length_elements + + ] = nextids [ i ] ;
// add new edges
nextids_count = 0 ;
for ( int i = elements ; i < elements + length_elements ; i + + )
for ( int j = 0 ; j < rank ; j + + ) {
edges [ i * rank + j ] = multiply_generator ( j , ids [ i ] , simple_roots , root_mapping , rank , positive ) ;
if ( ! ( ids [ i ] & BIT ( simple_roots [ j ] ) ) ) // the new element is longer then the old one
nextids [ nextids_count + + ] = edges [ i * rank + j ] ;
}
elements + = length_elements ;
}
// translate the ids to list positions (i.e. local continuous ids)
for ( int i = 0 ; i < order ; i + + ) {
lookup [ i ] . id = ids [ i ] ;
lookup [ i ] . position = i ;
}
qsort ( lookup , order , sizeof ( weylid_lookup_t ) , compare_weylid_lookup ) ;
for ( int i = 0 ; i < order ; i + + ) {
group [ i ] . id = ids [ i ] ;
for ( int j = 0 ; j < rank ; j + + )
group [ i ] . left [ j ] = lookup_id ( edges [ i * rank + j ] , lookup , order ) ;
}
free ( cartan_matrix ) ;
free ( root_vectors ) ;
free ( vector ) ;
free ( root_mapping ) ;
free ( simple_roots ) ;
free ( ids ) ;
free ( edges ) ;
free ( nextids ) ;
free ( lookup ) ;
}