#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "include.h"
#include "group.h"


/*  This file contains subroutines associated with the
    group structure of the lattice symmetries. 
    The first multiplets must be:
      the trivial subgroup
      even transverse parity (full)
      odd transverse parity  (full)
      singlet representation of the full group
      various other representations of the full group.  
*/

/*  subgroup and multiplet return a pointer to constant objects.
    Thus, other code can safely use pointers to specify subgroups
    and multiplets.
    Other code must not write onto these objects.   */


#if HINDEX == 1

/******************************************************************

     Multiplets are as follows:
     Charge conjugation: C = P_1 O, where O is orientation reversal
     In the following P_1 is parity and E is the identity.

   group   |   Multiplet                           
   element | 0  1  2  3   
  ---------|------------
   E       | 1  1  1  1  
   P_1     |    1 -1  1                      
  ----------------------

*******************************************************************/

const int subgroup[]={0,1,1,1};  /* subgroup for each multiplet */

/* 
    list of all lattice symmetry groups
    must have TRIVIAL_GROUP then FULL_GROUP 
*/

const latlist latgroups[]={{{0},1}, {{0,1},2}}; 

/*
   weight factors for multiplets
*/

const int multiplet[NMULTIPLETS][ORDER_GROUP]={{1},{1,1},{1,-1},{1,1}};

#elif HINDEX == 2


/*******************************************************************

     Multiplets are as follows:
     Charge conjugation: C = R^2 O, where O is orientation reversal
     In the following P_1 and P_2 are parity, R is a rotation,
     and E is the identity.
     Note that:  P_2 = P_1 R^2 and R^2 = P_1 P_2.

  group   |                     Multiplet                           
  element | 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 
 ---------|---------------------------------------------------------
  E       | 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
  P_1     |          1 -1  1 -1  1 -1  1 -1                         
  P_2     |          1 -1  1 -1 -1  1  1 -1              1 -1       
  R^2     |    1 -1  1  1  1  1 -1 -1  1  1  1 -1  1 -1             
          |                                                         
  P_1 R   |          1 -1 -1  1              1  1 -1 -1        1 -1 
  R       |          1  1 -1 -1                                     
  R^3     |          1  1 -1 -1                                     
  P_1 R^3 |          1 -1 -1  1              1 -1 -1  1             
 -------------------------------------------------------------------
  J_z^P_1 |         0+ 0- 2+ 2- 1+ 1-

  The multiplets have the following use:

  9,  8    P_perp = (c,0),   P_2 = +1
  7, 10    P_perp = (c,0),   P_2 = -1
  11,12    P_perp = (c,c), P_1 R = +1
  14,13    P_perp = (c,c), P_1 R = -1
  15       n = (c,0)         P_2 = +1        
  16       n = (c,0)         P_2 = -1
  17       n = (c,c)       P_1 R = +1 
  18       n = (c,c)       P_1 R = -1

********************************************************************/

const int subgroup[NMULTIPLETS]={                /* subgroup for each multiplet */
  TRIVIAL_GROUP,2,2,1,1,1,1,
  3,3,3,3, 4,4,4,4, 5,5, 6,6
};

/* 
    list of all lattice symmetry groups
    must have TRIVIAL_GROUP then FULL_GROUP first
*/

const latlist latgroups[]={
  {{0},1}, {{0,1,2,3,4,5,6,7},8},
  {{0,3},2}, {{0,1,2,3},4}, {{0,3,4,7},4},
  {{0,2},2}, {{0,4},2}
};

/*
   weight factors for multiplets
*/

const int multiplet[NMULTIPLETS][ORDER_GROUP]={
  {1}, {1,1}, {1,-1},

  {1,1,1,1, 1,1,1,1}, {1,-1,-1,1, -1,1,1,-1},   
  {1,1,1,1, -1,-1,-1,-1}, {1,-1,-1,1, 1,-1,-1,1},
  {1,1,-1,-1}, {1,-1,1,-1}, {1,1,1,1}, {1,-1,-1,1},

  {1,1,1,1}, {1,-1,1,-1}, {1,1,-1,-1}, {1,-1,-1,1},

  {1,1}, {1,-1}, {1,1}, {1,-1}
};

#else

#endif

/* permutations of HINDEX objects */

#if HINDEX==1
const int permutations[1]={0}; 
#elif HINDEX==2
const int permutations[4]={0,1,1,0}; 
#else
#endif

/*
    The function latorder permutes and flips the various winding number
    indices into a canonical order based on the lattice symmetry group, 
    returning the associated lattice symmetry group element
    needed to undo the transformation.  
*/

latype latorder(int *hlist){
  int i,h1[HINDEX],h2[HINDEX],h3[HINDEX];
  const int *perm;
  latype p,pbest=0;
#if 0  /* debug print */
  printf("        in (%i,%i)",hlist[0],hlist[1]);
#endif
  for(i=0; i<HINDEX; i++)h1[i]=hlist[i];
  for(p=1; p<ORDER_GROUP; p++){
    for(i=0; i<HINDEX; i++)h2[i]=((p&(1<<i))?-h1[i]:h1[i]);
    perm=permutations+(p>>HINDEX)*HINDEX;
    for(i=0; i<HINDEX; i++)h3[perm[i]]=h2[i];
#if 0 /* debug print */
    printf(" (%i,%i)",h3[0],h3[1]);
#endif
    for(i=0; i<HINDEX && h3[i]==hlist[i]; i++);
    if(i<HINDEX&&h3[i]>hlist[i]){                   
      for(i=0; i<HINDEX; i++)hlist[i]=h3[i];
      pbest=p;
    }  
  } 
#if 0  /* debug print */
  printf(" out (%i,%i) p=%i\n",hlist[0],hlist[1],pbest);
#endif
  return pbest;
}

/* Check to see if an integer vector is invariant
   under a given subgroup of the lattice rotation group.
   it returns a boolean value. */

int isinvariantvector(int *n,latlist *sub){
  int i,j,n2[HINDEX];
  const int *perm;
  for(i=0; i<sub->length; i++){
    perm=permutations+(sub->e[i]>>HINDEX)*HINDEX;
    for(j=0; j<HINDEX; j++)
      n2[perm[j]]=n[j];
    for(j=0; j<HINDEX; j++)
      n2[j]*=1-2*((sub->e[i]>>j)&1);
#if 0 /* debug print */
    printf("  element %i n=(%i,%i) n2=(%i,%i)\n",
	   sub->e[i],n[0],n[1],n2[0],n2[1]);
#endif
    for(j=0;j<HINDEX;j++)
      if(n[j]!=n2[j])return 0;
  }
  return 1;
}

/*This checks that the n is invariant.*/

int fisinvariantvector(element *n, latlist *sub){
  int i,j,k;
  const int *perm;
  element n2[HINDEX];
  for(i=0; i<sub->length; i++){
    perm=permutations+(sub->e[i]>>HINDEX)*HINDEX;
    for(j=0; j<HINDEX; j++)
      n2[perm[j]]=n[j];
    for(j=0; j<HINDEX; j++){
      k=(1-2*((sub->e[i]>>j)&1));
      n2[j]*=k;
    }
#if 1 /* debug print */
    printf("  element %i n=(%g,%g) n2=(%g,%g)\n",
	   sub->e[i],n[0],n[1],n2[0],n2[1]);
#endif
    for(j=0; j<HINDEX; j++)
    if(fabs(n[j]-n2[j]) > 2*ELEMENT_EPSILON)return 0;
  }
  return 1;
}