// Compute number of zero circulant determinant from {+1,-1} version 3; 
// C-language GNU C99 on Apple Powerbook G4, Fred Lunnon, Maynooth 04/10/08 
// Compile, run, save via:   cd /home/fred/macwall; 
//   gcc -c circulant.c; gcc -o circulant circulant.o -lm; 
//   time ./circulant 0 24; cp circulant.c circulant3.c; 

/* The program enumerates singular order n circulants with {-1,+1} components, 
where n < 64. No advantage is taken of known results, such as circ_count(n) = 
    2 for n odd prime; 2 n_C_n/2 for n twice odd prime; 
    -\sum_{d | n, d < n} \mu(n/d) \sum_{k = 0..d} (d_C_k)^(n/d) 
    for n odd square-free (Alekseyev, Jovovic). 

A CLI is provided accepting a list of low, high bounds of intervals of n, 
such as <circulant 0 24 32 32>. 

As S ranges over binary n-sequences, for i zero and each i dividing n, 
the algorithm tests whether  \sum_j  S_j w^{i j}  =  0  for 0 <= j < n; 
where w is a primitive n-th root of unity, and S_j = (+/-) 1. 
Time cost is order (n log n) <long long> (64-bit) additions per S 
for those that fail (the majority). 

Isomorphs of sequences under cycling, reversal, complementation are rejected, 
at a cost of order (4 n) bitwise logical ops per canonical S, ~1/(4 n) of cases. 
Only canonical S (numerically smaller than any isomorphs) are tested; 
such S carry a weight <syms> equal to the number of their isomorphs. 

Further, the outer loop is curtailed to avoid generating some isomorphs: 
S is forced to start 00... and finish ...1, while 00111... is forbidden; 
the number of strings generated is reduced thereby from 2^n to (7/64)2^n. 
Special cases S = 00...0 and 0101...01 are treated separately. 

Time for n = 0--24 on Apple Mac G4 Powerboook 
    8.5 = 0.5 (gen) + 6.5 (isom) + 1.5 (det) secs; for n = 32, 20 min; 
Time order (n log n) 2^n, storage order n log n. 
Use of <long long> to manipulate sequences limits n < 64. 
No warning or detection of integer overflow is implemented! 

The test might be modified to breakdown by rank, which might be of interest; 
also to enumerate singular {0,1} circulants instead, tho' may need to refrain 
from extracting complement symmetry! 

FFT might be used to evaluate these determinants for n = 2^k 
(also, with increasing difficulty, for other highly composite n). 

When a non-canonical S is rejected, perhaps advance over the block of S 
which will be rejected at the same point: e.g.  001011100... -> 001001001... 
The cycle loop should shift up/left to make this effective. */ 

#include <math.h> // sin(), cos(), fabs() 

// Globals set by circ_count(), accessed by isom_det() 
#define nmax 63 
int n, nsig; // circulant maximum, order, number of factors 
long long int pmmax, pmint; 
double eps = 1.0E-12; // linear factor approx. zero (using 17 sig figs) 
double pi; 
double rept [nmax][nmax], impt[nmax][nmax]; // cos and sin tables 

// Reject isomorphs and test singular circulant for row = pmint, order = n; 
//   [actually slightly faster if reverse rejection omitted!] 
int isom_det () {int i,j; 
  long long int pmrot, pmrev, pmcom; 
  //long long int pmmax; double eps; double rept [n][n], impt[n][n]; 
  //int n, nsig; // globals 
  int syms; double facr, faci; 

  // Reject unless canonical; count isomorphs 
  pmrot = pmint; syms = 0; 
  for (j = 0; j < n && pmrot >= pmint; j += 1) /*do*/ {// cycle 
    if (pmrot == pmint) /*then*/ syms += 1; /*fi*/ 
    pmrot = ((pmrot >> 1) ^ (pmrot << (n-1))) & pmmax; // (circ. shift) 
    } /*od*/ 
  if (pmrot < pmint) /*then*/ return(0); /*fi*/ 
  pmcom = (-1-pmint) & pmmax; // complement 
  for (j = 0; j < n && pmcom >= pmint; j += 1) /*do*/ {// cycle  
    if (pmcom == pmint) /*then*/ syms += 1; /*fi*/ 
    pmcom = ((pmcom >> 1) ^ (pmcom << (n-1))) & pmmax; 
    } /*od*/ 
  if (pmcom < pmint) /*then*/ return(0); /*fi*/ 

  pmrev = 0; // (pmrot == pmint) reverse 
  for (j = 0; j < n; j += 1) /*do*/ {
    pmrev = (pmrev << 1) | (pmrot & 1); 
    pmrot = (pmrot >> 1);} /*od*/ 

  for (j = 0; j < n && pmrev >= pmint; j += 1) /*do*/ {// cycle 
    if (pmrev == pmint) /*then*/ syms += 1; /*fi*/ 
    pmrev = ((pmrev >> 1) ^ (pmrev << (n-1))) & pmmax; // (circ. shift) 
    } /*od*/ 
  if (pmrev < pmint) /*then*/ return(0); /*fi*/ 
  pmcom = (-1-pmrev) & pmmax; // complement 
  for (j = 0; j < n && pmcom >= pmint; j += 1) /*do*/ {// cycle  
    if (pmcom == pmint) /*then*/ syms += 1; /*fi*/ 
    pmcom = ((pmcom >> 1) ^ (pmcom << (n-1))) & pmmax; 
    } /*od*/ 
  if (pmcom < pmint) /*then*/ return(0); /*fi*/ 

  // Detect zero linear factor of singular circulant 
  // [What is the optimal order in which to test these?] 
  for (i = 0; i < nsig; i += 1) /*do*/ {
    facr = 0.0; faci = 0.0; pmrot = pmint; 
    for (j = 0; j < n; j += 1) /*do*/ {
      if (pmrot & 1) /*then*/ 
        {facr = facr + rept[i][j]; 
        faci = faci + impt[i][j];} 
      else 
        {facr = facr - rept[i][j]; 
        faci = faci - impt[i][j];} 
      pmrot = pmrot >> 1; // shift down 
      } /*od*/ 
    if (fabs(facr) < eps && fabs(faci) < eps) /*then*/ 
      break; /*fi*/} /*od*/ // beware: abs() fails to cast! 

  return(/*if*/ (i < nsig) ?/*then*/ 4*n/syms :/*else*/ 0 /*fi*/);} 

// Count singular circulant determinants of order n over {-1,+1}, n < 64 
long long int circ_count (n0) 
  int n0; {
  //long long int pmmax, pmint; double eps, pi2, rept [n][n], impt[n][n]; 
  //int n, nsig; // globals 
  n = n0; 
  // Set up cos and sin tables for roots of unity in isom_det() 
  int i,j,k = 0; pi = 4.0*atan(1.0); // cannot set globally! 
  for (i = 0; i < n; i += 1) /*do*/ 
    if (i == 0 || n%i == 0) /*then*/ {for (j = 0; j < n; j += 1) /*do*/ 
      {rept[k][j] = cos(i*j*2*pi/n); impt[k][j] = sin(i*j*2*pi/n);} /*od*/ 
      k += 1;} /*fi od*/ 
  nsig = k; // record factor classes = divisors of n 
  long long int cter = 0; 
  pmmax = 1; pmmax = (pmmax << n) - 1; 
    // 11...1 for isom_det(); cast 64-bit shift! 

  // Extract special cases 
  pmint = 0; cter = cter + isom_det(); // 00...0 
  for (j = 0; j < n/2; j += 1) /*do*/ pmint = (pmint << 2) | 1; /*od*/ 
  cter = cter + isom_det(); // 0101..01 

  // Main loop: force first bits = 00, next != 111, last = 1 avoiding isomorphs; 
  long long int pmfin = 7; pmfin = (pmfin << (n-5)) - 3; // cast 64-bit shift! 
  if (n < 5) /*then*/ pmfin = (pmmax >> 2) - 2; /*fi*/ // small cases 
  pmint = -1; // beware: long long fails as for-loop control! 
  while (pmint <= pmfin) /*do*/ {pmint += 2; 
    cter = cter + isom_det();} /*od*/ 
  return(cter);} 

int main (argc, argv) 
  int argc; char **argv; {
  int i, na, nb, n; long long int ccn; 
  printf("   n            count \n"); 
  i = 0; while (i < argc-1) /* do */ {
    // input start and finish from CLI, disable parameters 
    na = atoi(argv[++i]); argv[i] = argv[0]; 
    nb = atoi(argv[++i]); argv[i] = argv[0]; 
    for (n = na; n <= nb; n += 1) /*do*/ {
      ccn = circ_count (n); 
      printf("%4d %16lld \n", n, ccn);} /*od*/} /*od*/ 
  return;} 

// 0, 0, 4, 2, 8, 2, 40, 2, 128, 62, 
// 504, 2, 1768, 2, 6864, 2738, 24192, 2, 107128, 2, 
// 331096, 109334, 1410864, 2, 5880544, 206282, 20801200, 5417630, 73508696, 2, 
// 345334744, 2, 1150681600, 278770214, 4667212440, 133401818, 19355912632, 2, 70690527600, 

/* 
   n        count 
   0            0 
   1            0 
   2            4 
   3            2 
   4            8 
   5            2 
   6           40 
   7            2 
   8          128 
   9           62 
  10          504 
  11            2 
  12         1768 
  13            2 
  14         6864 
  15         2738 
  16        24192 
  17            2 
  18       107128 
  19            2 
  20       331096 
  21       109334 
  22      1410864 
  23            2 
  24      5880544 
  25       206282 
  26     20801200 
  27      5417630 
  28     73508696 
  29            2 
  30    345334744 
  31            2 
  32   1150681600 
  33    278770214 
  34   4667212440 
  35    133401818 
  36  19355912632 
  37            2 
  38  70690527600 
  39  15151534406
*/ 
// 1027 mins for n = 33--36 version 2; approx 18 hrs for 0--36 (now 12?) 
// 37--38 via Miller & Alekseyev theorems; 39 Alekseyev unchecked?