/*
  Space Group Info's (c) 1994-96 Ralf W. Grosse-Kunstleve
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


#define SGCLIB_C__
#include "sginfo.h"


static const char *Err_Ill_SMx_in_List =
  "Error: Illegal SeitzMx in list";


void SetSgError(const char *msg)
{
  if (SgError == NULL) SgError = msg;
}


int iModPositive(int ix, int iy)
{
  if (iy > 0)
  {
    ix %= iy;
    if (ix < 0) ix += iy;
  }

  return ix;
}


static void SwapRTMx(T_RTMx *Mx_a, T_RTMx *Mx_b)
{
  int     i;
  T_RTMx  BufMx;


  for (i = 0; i < 12; i++) BufMx.a[i] = Mx_a->a[i];
  for (i = 0; i < 12; i++) Mx_a->a[i] = Mx_b->a[i];
  for (i = 0; i < 12; i++) Mx_b->a[i] = BufMx.a[i];
}


static void CopyRotMxInfo(T_RotMxInfo *target, const T_RotMxInfo *source)
{
  memcpy(target, source, sizeof (*target));
}


static void SwapRotMxInfo(T_RotMxInfo *RMx_a, T_RotMxInfo *RMx_b)
{
  T_RotMxInfo  Buffer;

  memcpy(&Buffer, RMx_a,   sizeof (Buffer));
  memcpy(RMx_a,   RMx_b,   sizeof (Buffer));
  memcpy(RMx_b,   &Buffer, sizeof (Buffer));
}


int traceRotMx(const int *RotMx)
{
  return RotMx[0] + RotMx[4] + RotMx[8];
}


int deterRotMx(const int *RotMx)
{
  int     det;

  det =  RotMx[0] * (RotMx[4] * RotMx[8] - RotMx[5] * RotMx[7]);
  det -= RotMx[1] * (RotMx[3] * RotMx[8] - RotMx[5] * RotMx[6]);
  det += RotMx[2] * (RotMx[3] * RotMx[7] - RotMx[4] * RotMx[6]);

  return det;
}


void RotMx_t_Vector(int *R_t_V, const int *RotMx, const int *Vector, int FacTr)
{
  const int  *vec;


  if (FacTr > 0)
  {
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec;
    *R_t_V   %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
     R_t_V++;
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec;
    *R_t_V   %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
     R_t_V++;
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V   += *RotMx   * *vec;
    *R_t_V   %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
  }
  else
  {
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V++ += *RotMx++ * *vec;
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V++ += *RotMx++ * *vec;
                            vec = Vector;
    *R_t_V   =  *RotMx++ * *vec++;
    *R_t_V   += *RotMx++ * *vec++;
    *R_t_V   += *RotMx   * *vec;
  }
}


void RotMxMultiply(int *rmxab, const int *rmxa, const int *rmxb)
{
  const int  *a, *b;

  /* no loops to be as fast as posslible */

  a = rmxa;
  b = rmxb;
  *rmxab  = *a++ * *b; b += 3; /* r11 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r12 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r13 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a++ * *b; b = rmxb;
   rmxab++;

  rmxa = a;
  *rmxab  = *a++ * *b; b += 3; /* r21 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r22 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r23 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a++ * *b; b = rmxb;
   rmxab++;

  rmxa = a;
  *rmxab  = *a++ * *b; b += 3; /* r31 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r32 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b; b -= 5;
   rmxab++;

  a = rmxa;
  *rmxab  = *a++ * *b; b += 3; /* r33 */
  *rmxab += *a++ * *b; b += 3;
  *rmxab += *a   * *b;
}


void RotateRotMx(int *RotMx, const int *RMx, const int *InvRMx)
{
  int  BufMx[9];


  RotMxMultiply(BufMx, RotMx, InvRMx);
  RotMxMultiply(RotMx, RMx,   BufMx);
}


void SeitzMxMultiply(T_RTMx *smxab, const T_RTMx *smxa, const T_RTMx *smxb)
{
  const int  *ar, *a, *b, *bt;
  int        *ab;

  /* no loops to be as fast as posslible */

  ar = smxa->a;
  a  = smxa->a;
  b  = smxb->a;
  ab = smxab->a;

  *ab  = *a++ * *b; b += 3; /* r11 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r12 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r13 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b; b = smxb->a;
   ab++;

  ar = a;
  *ab  = *a++ * *b; b += 3; /* r21 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r22 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r23 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b; b = smxb->a;
   ab++;

  ar = a;
  *ab  = *a++ * *b; b += 3; /* r31 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r32 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r33 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b++; bt = b;
   ab++;

  ar = smxa->a;
  *ab  = *ar++ * *b++; /* t1 */
  *ab += *ar++ * *b++;
  *ab += *ar++ * *b; b = bt;
  *ab += *a++;
  *ab %= STBF; if (*ab < 0) *ab += STBF;
   ab++;

  *ab  = *ar++ * *b++; /* t2 */
  *ab += *ar++ * *b++;
  *ab += *ar++ * *b; b = bt;
  *ab += *a++;
  *ab %= STBF; if (*ab < 0) *ab += STBF;
   ab++;

  *ab  = *ar++ * *b++; /* t3 */
  *ab += *ar++ * *b++;
  *ab += *ar   * *b;
  *ab += *a;
  *ab %= STBF; if (*ab < 0) *ab += STBF;
}


void RTMxMultiply(T_RTMx *rtmxab, const T_RTMx *rtmxa, const T_RTMx *rtmxb,
                  int FacAug, int FacTr)
{
  const int  *ar, *a, *b, *bt;
  int        *ab;

  /* no loops to be as fast as posslible */

  ar = rtmxa->a;
  a  = rtmxa->a;
  b  = rtmxb->a;
  ab = rtmxab->a;

  *ab  = *a++ * *b; b += 3; /* r11 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r12 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r13 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b; b = rtmxb->a;
   ab++;

  ar = a;
  *ab  = *a++ * *b; b += 3; /* r21 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r22 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r23 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b; b = rtmxb->a;
   ab++;

  ar = a;
  *ab  = *a++ * *b; b += 3; /* r31 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r32 */
  *ab += *a++ * *b; b += 3;
  *ab += *a   * *b; b -= 5;
   ab++;

  a = ar;
  *ab  = *a++ * *b; b += 3; /* r33 */
  *ab += *a++ * *b; b += 3;
  *ab += *a++ * *b++; bt = b;
   ab++;

  if (FacTr > 0)
  {
    ar = rtmxa->a;
    *ab  = *ar++ * *b++; /* t1 */
    *ab += *ar++ * *b++;
    *ab += *ar++ * *b; b = bt;
    *ab += *a++ * FacAug;
    *ab %= FacTr; if (*ab < 0) *ab += FacTr;
     ab++;

    *ab  = *ar++ * *b++; /* t2 */
    *ab += *ar++ * *b++;
    *ab += *ar++ * *b; b = bt;
    *ab += *a++ * FacAug;
    *ab %= FacTr; if (*ab < 0) *ab += FacTr;
     ab++;

    *ab  = *ar++ * *b++; /* t3 */
    *ab += *ar++ * *b++;
    *ab += *ar   * *b;
    *ab += *a   * FacAug;
    *ab %= FacTr; if (*ab < 0) *ab += FacTr;
  }
  else
  {
    ar = rtmxa->a;
    *ab  = *ar++ * *b++; /* t1 */
    *ab += *ar++ * *b++;
    *ab += *ar++ * *b; b = bt;
    *ab += *a++ * FacAug;
     ab++;

    *ab  = *ar++ * *b++; /* t2 */
    *ab += *ar++ * *b++;
    *ab += *ar++ * *b; b = bt;
    *ab += *a++ * FacAug;
     ab++;

    *ab  = *ar++ * *b++; /* t3 */
    *ab += *ar++ * *b++;
    *ab += *ar   * *b;
    *ab += *a   * FacAug;
  }
}


void InverseRotMx(const int *RotMx, int *InvRotMx)
{
  InvRotMx[0] =   RotMx[4] * RotMx[8] - RotMx[5] * RotMx[7];
  InvRotMx[1] = - RotMx[1] * RotMx[8] + RotMx[2] * RotMx[7];
  InvRotMx[2] =   RotMx[1] * RotMx[5] - RotMx[2] * RotMx[4];
  InvRotMx[3] = - RotMx[3] * RotMx[8] + RotMx[5] * RotMx[6];
  InvRotMx[4] =   RotMx[0] * RotMx[8] - RotMx[2] * RotMx[6];
  InvRotMx[5] = - RotMx[0] * RotMx[5] + RotMx[2] * RotMx[3];
  InvRotMx[6] =   RotMx[3] * RotMx[7] - RotMx[4] * RotMx[6];
  InvRotMx[7] = - RotMx[0] * RotMx[7] + RotMx[1] * RotMx[6];
  InvRotMx[8] =   RotMx[0] * RotMx[4] - RotMx[1] * RotMx[3];
}


void InverseRTMx(const T_RTMx *RTMx, T_RTMx *InvRTMx)
{
  int        *iR;
  const int  *T;


  iR = InvRTMx->s.R;

  InverseRotMx(RTMx->s.R, iR);

  T = RTMx->s.T;

  InvRTMx->s.T[0] = - iR[0] * T[0] - iR[1] * T[1] - iR[2] * T[2];
  InvRTMx->s.T[1] = - iR[3] * T[0] - iR[4] * T[1] - iR[5] * T[2];
  InvRTMx->s.T[2] = - iR[6] * T[0] - iR[7] * T[1] - iR[8] * T[2];
}


int IsSMxTransl0(const T_LatticeInfo *LatticeInfo, const int *SeitzMxT)
{
  int        iTrV, nTrV, t;
  const int  *TrV;


  nTrV = LatticeInfo->nTrVector;
   TrV = LatticeInfo->TrVector;

  for (iTrV = 0; iTrV < nTrV; iTrV++)
  {
        t =     (SeitzMxT[0] + TrV[0]) % STBF;
    if (t == 0) {
          t =   (SeitzMxT[1] + TrV[1]) % STBF;
      if (t == 0) {
            t = (SeitzMxT[2] + TrV[2]) % STBF;
        if (t == 0)
          return 1;
    }}

    TrV += 3;
  }

  return 0;
}


static int IsSpecialSeitzMx(T_SgInfo *SgInfo, const T_RTMx *SMx, int ExpandLT)
{
  int                  i, special, smx11;
  const T_LatticeInfo  *ExpLT;


  /* check rotation part for identity or inversion operation */

         smx11 = SMx->s.R[0];
  if (   smx11 !=  1
      && smx11 != -1) return 0;

  for (i = 1; i < 9; i++)
  {
    if (i % 4) {
      if (SMx->s.R[i] !=     0) return 0; }
    else {
      if (SMx->s.R[i] != smx11) return 0; }
  }

  if (smx11 == 1) special = SpecialSMx_Identity;
  else            special = SpecialSMx_Inversion;

  /* rotation part is identity or inversion
     check translation part now
   */

  if (IsSMxTransl0(SgInfo->LatticeInfo, SMx->s.T) == 1)
    return (special | SpecialSMx_Transl0);

  if (ExpandLT && (smx11 == 1 || SgInfo->Centric))
  {
    /* try to expand lattice type */

    ExpLT = NULL;

    switch (SgInfo->LatticeInfo->Code)
    {
      case 'P':
        if (IsSMxTransl0(LI_A, SMx->s.T) == 1)
        { ExpLT = LI_A; break; }
        if (IsSMxTransl0(LI_B, SMx->s.T) == 1)
        { ExpLT = LI_B; break; }
        if (IsSMxTransl0(LI_C, SMx->s.T) == 1)
        { ExpLT = LI_C; break; }
        if (IsSMxTransl0(LI_I, SMx->s.T) == 1)
        { ExpLT = LI_I; break; }
        if (IsSMxTransl0(LI_R, SMx->s.T) == 1)
        { ExpLT = LI_R; break; }
        if (IsSMxTransl0(LI_S, SMx->s.T) == 1)
        { ExpLT = LI_S; break; }
        if (IsSMxTransl0(LI_T, SMx->s.T) == 1)
        { ExpLT = LI_T; break; }
      case 'A':
      case 'B':
      case 'C':
        if (IsSMxTransl0(LI_F, SMx->s.T) == 1)
          ExpLT = LI_F;
    }

    if (ExpLT != NULL)
    {
      SgInfo->LatticeInfo = ExpLT;
      SgInfo->StatusLatticeTr = -1;
      return (special | SpecialSMx_Transl0);
    }
  }

  return special;
}


int CompareSeitzMx(const T_LatticeInfo *LatticeInfo,
                   const T_RTMx *SeitzMxA, const T_RTMx *SeitzMxB)
{
  int        i, iTrV, nTrV, t;
  const int  *TrV;


  /* compare rotation part */

  for (i = 0; i < 9; i++)
    if (SeitzMxA->s.R[i] != SeitzMxB->s.R[i]) return 1;

  /* rotation part is same
     check translation
   */

  nTrV = LatticeInfo->nTrVector;
   TrV = LatticeInfo->TrVector;

  for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
  {
        t =     (SeitzMxA->s.T[0] + TrV[0]) % STBF;
    if (t ==     SeitzMxB->s.T[0]) {
          t =   (SeitzMxA->s.T[1] + TrV[1]) % STBF;
      if (t ==   SeitzMxB->s.T[1]) {
            t = (SeitzMxA->s.T[2] + TrV[2]) % STBF;
        if (t == SeitzMxB->s.T[2])
          return 0;
    }}
  }

  return -1;
}


int GetRotMxOrder(const int *RotMx)
{
  int deter = deterRotMx(RotMx);

  if (deter == -1 || deter == 1)
  {
    switch (traceRotMx(RotMx))
    {
      case -3:                  return -1;
      case -2:                  return -6;
      case -1: if (deter == -1) return -4;
               else             return  2;
      case  0: if (deter == -1) return -3;
               else             return  3;
      case  1: if (deter == -1) return -2;
               else             return  4;
      case  2:                  return  6;
      case  3:                  return  1;
    }
  }

  return 0;
}


static int nNextBasis_of_DirCode(const int DirCode,
                                 const int **RMx, const int **InvRMx)
{
  switch (DirCode)
  {
    case  '.': *RMx =            *InvRMx = NULL;      return 1;
    case  '=':
    case  '"':
    case '\'':
    case  '|':
    case '\\': *RMx = RMx_3_111; *InvRMx = RMx_3i111; return 3;
    case  '*': *RMx = RMx_4_001; *InvRMx = RMx_4i001; return 4;
    default:
      break;
  }

  SetSgError("Internal Error: Corrupt DirCode");
  return -1;
}


int GetRotMxInfo(const int *RotMx, T_RotMxInfo *RotMxInfo)
{
  int        i;
  int        nNextBasis, iNextBasis;
  int        nLoopInv, iLoopInv;
  int        Order, AbsOrder;
  int        RMxCopy[9], MatchMx[9], InvMatchMx[9], REV[3];
  int        *mmx;
  const int  *NBRMx, *InvNBRMx;

  const T_TabXtalRotMx  *txrmx;


  Order = GetRotMxOrder(RotMx);

  if (RotMxInfo)
      RotMxInfo->Order = Order;

  if (Order)
  {
    AbsOrder = abs(Order);

    if (Order > 0)
      for (i = 0; i < 9; i++) RMxCopy[i] =  RotMx[i];
    else
      for (i = 0; i < 9; i++) RMxCopy[i] = -RotMx[i];

    for (txrmx = TabXtalRotMx; txrmx->Order; txrmx++)
      if (txrmx->Order == AbsOrder) break;

    while (txrmx->Order == AbsOrder)
    {
      nNextBasis = nNextBasis_of_DirCode(txrmx->DirCode, &NBRMx, &InvNBRMx);

      if (nNextBasis < 0)
        return 0;

      if (AbsOrder > 2) nLoopInv = 2;
      else              nLoopInv = 1;

      for (i = 0; i < 9; i++) MatchMx[i] = txrmx->RMx[i];

      for (iNextBasis = 0; iNextBasis < nNextBasis; iNextBasis++)
      {
        if (iNextBasis)
          RotateRotMx(MatchMx, NBRMx, InvNBRMx);

        mmx = MatchMx;

        for (iLoopInv = 0; iLoopInv < nLoopInv; iLoopInv++)
        {
          if (iLoopInv)
          {
            InverseRotMx(MatchMx, InvMatchMx);
            mmx = InvMatchMx;
          }

          for (i = 0; i < 9; i++)
            if (mmx[i] != RMxCopy[i]) break;

          if (i == 9) /* matrix has been found */
          {
            if (RotMxInfo)
            {
              RotMxInfo->Inverse = iLoopInv;

              if (nNextBasis == 3)
              {
                switch(iNextBasis)
                {
                  case 0: RotMxInfo->RefAxis = 'z'; break;
                  case 1: RotMxInfo->RefAxis = 'x'; break;
                  case 2: RotMxInfo->RefAxis = 'y'; break;
                }
              }
              else
                RotMxInfo->RefAxis = 'o';

              RotMxInfo->DirCode = txrmx->DirCode;

              for (i = 0; i < 3; i++)
                RotMxInfo->EigenVector[i] = txrmx->EigenVector[i];

              for (; iNextBasis--;)
              {
                RotMx_t_Vector(REV, NBRMx, RotMxInfo->EigenVector, 0);

                if (iNextBasis-- == 0)
                {
                  for (i = 0; i < 3; i++)
                    RotMxInfo->EigenVector[i] = REV[i];

                  break;
                }

                RotMx_t_Vector(RotMxInfo->EigenVector, NBRMx, REV, 0);
              }
            }

            return Order;
          }
        }
      }

      txrmx++;
    }
  }

  return 0;
}


const T_RotMxInfo *ListOrBufRotMxInfo(const T_SgInfo *SgInfo, int iList,
                                      T_RotMxInfo *BufRotMxInfo)
{
  T_RotMxInfo  *RMxI;


      RMxI = SgInfo->ListRotMxInfo;
  if (RMxI)
      RMxI += iList;
  else
  {
    RMxI = BufRotMxInfo;

    if (GetRotMxInfo(SgInfo->ListSeitzMx[iList].s.R, RMxI) == 0) {
      SetSgError(Err_Ill_SMx_in_List);
      return NULL;
    }
  }

  return RMxI;
}


static int CoreAdd2ListSeitzMx(T_SgInfo *SgInfo, const T_RTMx *NewSMx)
{
  int                  i, iList;
  T_RTMx               *lsmx;
  T_RotMxInfo          RotMxInfo;
  const T_LatticeInfo  *LI;

  static const char *Err_NonXtalOp =
    "Error: Generators produce non-crystallographic operation";


  if (SgInfo->GenOption) LI = SgInfo->LatticeInfo;
  else                   LI = LI_P;

  lsmx = SgInfo->ListSeitzMx;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
    if (CompareSeitzMx(LI, NewSMx, lsmx) == 0)
      return 0; /* matrix is not unique */

  if (GetRotMxInfo(NewSMx->s.R, &RotMxInfo) == 0) {
    SetSgError(Err_NonXtalOp);
    return -1;
  }

  if (SgInfo->nList >= SgInfo->MaxList)
  {
    if (SgInfo->nList >= 192)
      SetSgError(Err_NonXtalOp);
    else
      SetSgError("Internal Error: Allocated space for ListSeitzMx too small");

    return -1;
  }

  for (i = 0; i < 12; i++) lsmx->a[i] = NewSMx->a[i];

  if (SgInfo->ListRotMxInfo != NULL)
    CopyRotMxInfo(&SgInfo->ListRotMxInfo[SgInfo->nList], &RotMxInfo);

  SgInfo->nList++;

  return 1;
}


int Add2ListSeitzMx(T_SgInfo *SgInfo, const T_RTMx *NewSMx)
{
  int     i, special, Enter;
  int     iMult, jMult;
  T_RTMx  TrialSMx, *iMultSMx, *jMultSMx;

  static const char *Err_Ill_lattice_translation =
    "Error: Illegal lattice translation";


  if (SgInfo->nList == 0)
  {
    /* make sure identity matrix is first in list */

    InitSeitzMx(&TrialSMx, 1);

    if (CoreAdd2ListSeitzMx(SgInfo, &TrialSMx) < 0)
      return -1;;
  }

  for (i = 0; i < 9; i++)
        TrialSMx.s.R[i] = NewSMx->s.R[i];

  for (i = 0; i < 3; i++) {
        TrialSMx.s.T[i] = NewSMx->s.T[i] % STBF;
    if (TrialSMx.s.T[i] < 0)
        TrialSMx.s.T[i] += STBF;
  }

  iMult = SgInfo->nList;
  iMultSMx = &SgInfo->ListSeitzMx[iMult];

  jMult = 1;
  jMultSMx = &SgInfo->ListSeitzMx[1]; /* skip first = identity matrix */

  for (;;)
  {
    Enter = 1;

    special = IsSpecialSeitzMx(SgInfo, &TrialSMx, 1);

    if      (special & SpecialSMx_Identity)
    {
      if (! (special & SpecialSMx_Transl0)) {
        SetSgError(Err_Ill_lattice_translation);
        return -1;
      }

      if (SgInfo->GenOption)
        Enter = 0;
    }
    else if (special & SpecialSMx_Inversion)
    {
      if (! (special & SpecialSMx_Transl0))
      {
        if (SgInfo->Centric) {
          SetSgError(Err_Ill_lattice_translation);
          return -1;
        }

        SgInfo->InversionOffOrigin = 1;
      }
      else
      {
        if (SgInfo->InversionOffOrigin)
          SgInfo->Centric = 1;

        SgInfo->InversionOffOrigin = 0;

        if (SgInfo->GenOption)
        {
          if (SgInfo->Centric == 0)
              SgInfo->Centric = -1;

          Enter = 0;
        }
        else
          SgInfo->Centric = 1;
      }
    }

    if (Enter && CoreAdd2ListSeitzMx(SgInfo, &TrialSMx) < 0)
      return -1;

    if (SgInfo->GenOption < 0)
      break;

    if (jMult > iMult)
    {
      iMult++;
      iMultSMx++;

      jMult = 1;
      jMultSMx = &SgInfo->ListSeitzMx[1]; /* skip first = identity matrix */
    }

    if (iMult == SgInfo->nList)
      break;

    SeitzMxMultiply(&TrialSMx, jMultSMx, iMultSMx);

    jMult++;
    jMultSMx++;
  }

  return 0;
}


int AddInversion2ListSeitzMx(T_SgInfo *SgInfo)
{
  T_RTMx  SeitzMx;

  InitSeitzMx(&SeitzMx, -1);
  return Add2ListSeitzMx(SgInfo, &SeitzMx);
}


int AddLatticeTr2ListSeitzMx(T_SgInfo *SgInfo,
                             const T_LatticeInfo *LatticeInfo)
{
  int        iTrV, nTrV;
  const int  *TrV;
  T_RTMx     SeitzMx;


  InitRotMx(SeitzMx.s.R, 1);

  nTrV = LatticeInfo->nTrVector;
   TrV = &LatticeInfo->TrVector[3]; /* skip first vector which is always 000 */

  for (iTrV = 1; iTrV < nTrV; iTrV++)
  {
    SeitzMx.s.T[0] = *TrV++;
    SeitzMx.s.T[1] = *TrV++;
    SeitzMx.s.T[2] = *TrV++;

    if (Add2ListSeitzMx(SgInfo, &SeitzMx) < 0)
      return -1;
  }

  if (SgInfo->GenOption)
    SgInfo->StatusLatticeTr = 0;
  else
    SgInfo->StatusLatticeTr = 1;

  return 0;
}


static int RemoveLatticeTr(T_SgInfo *SgInfo)
{
  int          iList, jList, i;
  T_RTMx       *smxi, *smxj, *lastsmx;
  T_RotMxInfo  *lrmxi, *lastrmxi;


  if (SgInfo->LatticeInfo->Code == 'P')
    return 0;

  if (SgInfo->StatusLatticeTr == -1)
  {
    if (AddLatticeTr2ListSeitzMx(SgInfo, SgInfo->LatticeInfo) < 0)
      return -1;
  }

  iList = 0;

  while (iList < SgInfo->nList)
  {
    smxi = &SgInfo->ListSeitzMx[iList];
    jList = iList + 1;

    while (jList < SgInfo->nList)
    {
      smxj = &SgInfo->ListSeitzMx[jList];

      if (CompareSeitzMx(SgInfo->LatticeInfo, smxi, smxj) == 0)
      {
        /* copy last element to this place */

        SgInfo->nList--;
        lastsmx = &SgInfo->ListSeitzMx[SgInfo->nList];
        for (i = 0; i < 12; i++) smxj->a[i] = lastsmx->a[i];

        if (SgInfo->ListRotMxInfo != NULL)
        {
          lrmxi =    &SgInfo->ListRotMxInfo[jList];
          lastrmxi = &SgInfo->ListRotMxInfo[SgInfo->nList];
          CopyRotMxInfo(lrmxi, lastrmxi);
        }
      }
      else
        jList++;
    }

    iList++;
  }

  SgInfo->StatusLatticeTr = 0;

  return 0;
}


static int RemoveInversion(T_SgInfo *SgInfo)
{
  int          i, iList, special, deter, Centric, InversionOffOrigin;
  T_RTMx       *lsmx, *smx, ProperSMx;
  T_RotMxInfo  *lrmxi, *rmxi;


  Centric = 0;
  InversionOffOrigin = 0;

  lsmx = SgInfo->ListSeitzMx;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
  {
    special = IsSpecialSeitzMx(SgInfo, lsmx, 0);

    if (special & SpecialSMx_Inversion)
    {
      if (special & SpecialSMx_Transl0)
        Centric = 1;
      else
        InversionOffOrigin = 1;

      break;
    }
  }

  if (InversionOffOrigin && Centric) {
    SetSgError("Internal Error: Corrupt SgInfo");
    return -1;
  }

  if (Centric == 0)
  {
    if (InversionOffOrigin) {
      SgInfo->Centric = 0;
      SgInfo->InversionOffOrigin = 1;
    }
    else
    {
      if (SgInfo->Centric) SgInfo->Centric = -1;
      SgInfo->InversionOffOrigin = 0;
    }
  }
  else
  {
    SgInfo->InversionOffOrigin = 0;

    lsmx  = SgInfo->ListSeitzMx;
    lrmxi = SgInfo->ListRotMxInfo;
    iList = 0;

    while (iList < SgInfo->nList)
    {
      deter = deterRotMx(lsmx->s.R);

      if (deter == -1 && SgInfo->Centric == -1)
      {
        for (i = 0; i < 9; i++)
              ProperSMx.s.R[i] = -lsmx->s.R[i];

        for (i = 0; i < 3; i++) {
              ProperSMx.s.T[i] = -lsmx->s.T[i] % STBF;
          if (ProperSMx.s.T[i] < 0)
              ProperSMx.s.T[i] += STBF;
        }

        smx = SgInfo->ListSeitzMx;

        for (i = 0; i < SgInfo->nList; i++, smx++)
          if (CompareSeitzMx(LI_P, &ProperSMx, smx) == 0) break;

        if (i == SgInfo->nList)
        {
          for (i = 0; i < 12; i++) lsmx->a[i] = ProperSMx.a[i];

          deter = deterRotMx(lsmx->s.R);

          if (deter != 1 || (lrmxi && GetRotMxInfo(lsmx->s.R, lrmxi) == 0)) {
            SetSgError(Err_Ill_SMx_in_List);
            return -1;
          }
        }
      }

      if      (deter == -1)
      {
        /* copy last element to this place */

            SgInfo->nList--;
        if (SgInfo->nList != iList)
        {
          smx = &SgInfo->ListSeitzMx[SgInfo->nList];
          for (i = 0; i < 12; i++) lsmx->a[i] = smx->a[i];

          if (lrmxi)
          {
            rmxi = &SgInfo->ListRotMxInfo[SgInfo->nList];
            CopyRotMxInfo(lrmxi, rmxi);
          }
        }
      }
      else if (deter ==  1)
      {
        lsmx++;
        if (lrmxi != NULL) lrmxi++;
        iList++;
      }
      else {
        SetSgError(Err_Ill_SMx_in_List);
        return -1;
      }
    }

    SgInfo->Centric = -1;
  }

  return 0;
}


int ApplyOriginShift(T_SgInfo *SgInfo)
{
  int     HaveOrSh, iList, i;
  T_RTMx  *lsmx, SMx;
  int     OrSh[3], lo[3];


  HaveOrSh = 0;

  for (i = 0; i < 3; i++) {
        OrSh[i] = SgInfo->OriginShift[i] * (STBF / 12);
    if (OrSh[i]) HaveOrSh = 1;
  }

  if (HaveOrSh == 0)
    return 0;

  lsmx = SgInfo->ListSeitzMx;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
  {
    RotMx_t_Vector(lo, lsmx->s.R, OrSh, STBF);

    for (i = 0; i < 3; i++)
      lsmx->s.T[i] = iModPositive(lsmx->s.T[i] - lo[i] + OrSh[i], STBF);
  }

  if (SgInfo->Centric == -1)
  {
    lsmx = &SMx;

    InitSeitzMx(lsmx, -1);

    RotMx_t_Vector(lo, lsmx->s.R, OrSh, STBF);

    for (i = 0; i < 3; i++)
      lsmx->s.T[i] = iModPositive(lsmx->s.T[i] - lo[i] + OrSh[i], STBF);

    if (CoreAdd2ListSeitzMx(SgInfo, lsmx) < 0)
      return -1;

    SgInfo->Centric = 0;
    SgInfo->InversionOffOrigin = 1;
  }

  return 1;
}


static void TidyTranslation(T_SgInfo *SgInfo)
{
  int        iList;
  int        iTrV, nTrV;
  const int  *TrV;
  T_RTMx     *lsmx;
  int        t1, t2, t3, *lt1, *lt2, *lt3, mint1, mint2, mint3;
  int        n0t, n0mint;


  nTrV = SgInfo->LatticeInfo->nTrVector;

  lsmx  = SgInfo->ListSeitzMx;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
  {
    mint1 = *(lt1 = &lsmx->s.T[0]);
    mint2 = *(lt2 = &lsmx->s.T[1]);
    mint3 = *(lt3 = &lsmx->s.T[2]);

    TrV = SgInfo->LatticeInfo->TrVector;

    for (iTrV = 0; iTrV < nTrV; iTrV++)
    {
      t1 = ((*lt1) + *TrV++) % STBF;
      t2 = ((*lt2) + *TrV++) % STBF;
      t3 = ((*lt3) + *TrV++) % STBF;

      n0t = (t1 == 0) + (t2 == 0) + (t3 == 0);
      n0mint = (mint1 == 0) + (mint2 == 0) + (mint3 == 0);

      if (    n0t > n0mint
          || (   n0t == n0mint
              && (   t1 + t2 + t3 <  mint1 + mint2 + mint3
                  || (   t1 + t2 + t3 == mint1 + mint2 + mint3
                      && (t1 < mint1 || (t1 == mint1 && t2 < mint2)))))) {
        mint1 = t1; mint2 = t2; mint3 = t3;
      }
    }

    *lt1 = mint1;
    *lt2 = mint2;
    *lt3 = mint3;
  }
}


static T_SgInfo *Pt_SgInfo_ListSortFunction = NULL;

static int SgInfoListSortFunction(const int *iList_a, const int *iList_b)
{
  int          val_a, val_b, n0_a, n0_b, i;
  T_RTMx       *smx_a, *smx_b;
  T_RotMxInfo  RotMxInfo_a, RotMxInfo_b, *rmxi_a, *rmxi_b;
  T_SgInfo     *SgInfo = Pt_SgInfo_ListSortFunction;


  if (SgError != NULL) return 0;

  if (SgInfo->ListRotMxInfo == NULL)
  {
    rmxi_a = &RotMxInfo_a;
    rmxi_b = &RotMxInfo_b;

    smx_a = &SgInfo->ListSeitzMx[*iList_a];
    smx_b = &SgInfo->ListSeitzMx[*iList_b];

    if (   GetRotMxInfo(smx_a->s.R, rmxi_a) == 0
        || GetRotMxInfo(smx_b->s.R, rmxi_b) == 0)
    {
      SetSgError(Err_Ill_SMx_in_List);
      return 0;
    }
  }
  else
  {
    rmxi_a = &SgInfo->ListRotMxInfo[*iList_a];
    rmxi_b = &SgInfo->ListRotMxInfo[*iList_b];
  }

  val_a = abs(rmxi_a->Order);
  val_b = abs(rmxi_b->Order);

  if (val_a == 1 && val_b != 1) return -1;
  if (val_a != 1 && val_b == 1) return  1;
  if (rmxi_a->Order == 1 && rmxi_b->Order != 1) return -1;
  if (rmxi_a->Order != 1 && rmxi_b->Order == 1) return  1;

  if (val_a != 1)
  {
    if (val_a > val_b) return -1;
    if (val_a < val_b) return  1;
    if (rmxi_a->Order > rmxi_b->Order) return -1;
    if (rmxi_a->Order < rmxi_b->Order) return  1;
  }

  n0_a = n0_b = 0;

  for (i = 0; i < 3; i++)
  {
    if (rmxi_a->EigenVector[i] == 0) n0_a++;
    if (rmxi_b->EigenVector[i] == 0) n0_b++;
  }
  if (n0_a > n0_b) return -1;
  if (n0_a < n0_b) return  1;

  val_a = val_b = 0;

  for (i = 0; i < 3; i++)
  {
    if (val_a < abs(rmxi_a->EigenVector[i]))
        val_a = abs(rmxi_a->EigenVector[i]);
    if (val_b < abs(rmxi_b->EigenVector[i]))
        val_b = abs(rmxi_b->EigenVector[i]);
  }
  if (val_a < val_b) return -1;
  if (val_a > val_b) return  1;

  val_a =  100 * abs(rmxi_a->EigenVector[2]);
  val_a +=  10 * abs(rmxi_a->EigenVector[0]);
  val_a +=       abs(rmxi_a->EigenVector[1]);
  val_b =  100 * abs(rmxi_b->EigenVector[2]);
  val_b +=  10 * abs(rmxi_b->EigenVector[0]);
  val_b +=       abs(rmxi_b->EigenVector[1]);

  if (n0_a < 2)
  {
    if (val_a < val_b) return -1;
    if (val_a > val_b) return  1;
  }
  else
  {
    if (val_a > val_b) return -1;
    if (val_a < val_b) return  1;
  }

  for (i = 0; i < 3; i++)
  {
    if (rmxi_a->EigenVector[i] > rmxi_b->EigenVector[i]) return -1;
    if (rmxi_a->EigenVector[i] < rmxi_b->EigenVector[i]) return  1;
  }

  if (rmxi_a->Inverse < rmxi_b->Inverse) return -1;
  if (rmxi_a->Inverse > rmxi_b->Inverse) return  1;

  smx_a = &SgInfo->ListSeitzMx[*iList_a];
  smx_b = &SgInfo->ListSeitzMx[*iList_b];

  for (i = 0; i < 3; i++)
  {
    if (smx_a->s.T[i] < smx_b->s.T[i]) return -1;
    if (smx_a->s.T[i] > smx_b->s.T[i]) return  1;
  }

  return 0;
}


static int PostSortSgInfoList(const T_SgInfo *SgInfo, int *List_iList)
{
  int                nList, iL_iL, jL_iL;
  T_RTMx             BufMx1, BufMx2, *smxab;
  const T_RTMx       *lsmx, *smxb;
  T_RotMxInfo        RotMxInfo;
  const T_RotMxInfo  *rmxi;
  int                nO_1, iO, save, i, i_;


  nList = SgInfo->nList;

  iL_iL = 0;

  while (iL_iL < nList)
  {
    lsmx = &SgInfo->ListSeitzMx[List_iList[iL_iL]];

        rmxi = ListOrBufRotMxInfo(SgInfo, List_iList[iL_iL], &RotMxInfo);
    if (rmxi == NULL)
      return -1;

    iL_iL++;

    iO = rmxi->Order;
    if (iO < 0 && iO % 2) iO *= 2;
    nO_1 = abs(iO) - 1;

    smxab = &BufMx2;
    smxb = lsmx;

    for (iO = 1; iO < nO_1; iO++)
    {
      SeitzMxMultiply(smxab, lsmx, smxb);

      for (jL_iL = iL_iL; jL_iL < nList; jL_iL++)
      {
        smxb = &SgInfo->ListSeitzMx[List_iList[jL_iL]];
        if (CompareSeitzMx(SgInfo->LatticeInfo, smxab, smxb) == 0)
          break;
      }

      if (jL_iL < nList)
      {
                              save = List_iList[jL_iL];

        for (i = i_ = jL_iL; i > iL_iL; i = i_)
          List_iList[i] = List_iList[--i_];

        List_iList[iL_iL++] = save;
      }

      smxb = smxab;
      if (iO % 2) smxab = &BufMx1;
      else        smxab = &BufMx2;
    }
  }

  return 0;
}


static void SortSgInfoList(T_SgInfo *SgInfo, int *List_iList)
{
  int          i, j, refi;
  int          nList;
  T_RTMx       *lsmx;
  T_RotMxInfo  *lrmxi;


  if (SgError != NULL) return;

  nList = SgInfo->nList;
  lsmx  = SgInfo->ListSeitzMx;
  lrmxi = SgInfo->ListRotMxInfo;
  Pt_SgInfo_ListSortFunction = SgInfo;

  for (i = 0; i < nList; i++) List_iList[i] = i;

  qsort((void *) List_iList, nList, sizeof (*List_iList),
        (int (*)(const void *, const void *)) SgInfoListSortFunction);

  Pt_SgInfo_ListSortFunction = NULL;
  if (SgError != NULL) return;

  if (PostSortSgInfoList(SgInfo, List_iList) != 0)
    return;

  for (i = 0; i < nList; i++)
  {
    j = List_iList[i];

    if (j != i)
    {
      for (refi = i + 1; refi < nList; refi++)
        if (List_iList[refi] == i) break;

      if (refi >= nList) {
        SetSgError("Internal Error: SortSgInfoList(): Corrupt List_iList");
        return;
      }

      SwapRTMx(&lsmx[i], &lsmx[j]);
      if (lrmxi != NULL)
        SwapRotMxInfo(&lrmxi[i], &lrmxi[j]);

      List_iList[refi] = j;
    }
  }
}


int FindSeitzMx(const T_SgInfo *SgInfo,
                int Order, int HonorSign, int RefAxis, int DirCode)
{
  int          iList;
  int          MatchingOrder;
  T_RTMx       *lsmx;
  T_RotMxInfo  *lrmxi, RotMxInfo;


  lsmx  = SgInfo->ListSeitzMx;
  lrmxi = SgInfo->ListRotMxInfo;

  if (lrmxi == NULL) lrmxi = &RotMxInfo;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
  {
    if (lrmxi == &RotMxInfo)
    {
      if (GetRotMxInfo(lsmx->s.R, lrmxi) == 0) {
        SetSgError(Err_Ill_SMx_in_List);
        return -1;
      }
    }

    if (HonorSign == 0)
      MatchingOrder = (abs(Order) == abs(lrmxi->Order));
    else
      MatchingOrder = (Order == lrmxi->Order);

    if (   MatchingOrder
        && lrmxi->Inverse == 0
        && (RefAxis == 0 || RefAxis == lrmxi->RefAxis)
        && (DirCode == 0 || DirCode == lrmxi->DirCode))
    {
      if (DirCode != '*') return iList;

      if (   lrmxi->EigenVector[0] == 1
          && lrmxi->EigenVector[1] == 1
          && lrmxi->EigenVector[2] == 1) return iList;
    }

    if (lrmxi != &RotMxInfo) lrmxi++;
  }

  return -1;
}


static int FindXtalSystem(T_SgInfo *SgInfo)
{
  int                iList, i;
  int                HonorSign = 0, CheckEnantiomorph;
  const T_RTMx       *lsmx;
  T_RotMxInfo        RotMxInfo;
  const T_RotMxInfo  *lrmxi;

  enum Axis { N1 = 0, N2, N3, N4, N6, EndOfAxis };
  int                         N_count[EndOfAxis];


  SgInfo->XtalSystem    = XS_Unknown;
  SgInfo->UniqueRefAxis = 0;
  SgInfo->UniqueDirCode = 0;
  SgInfo->ExtraInfo     = EI_Unknown;

  CheckEnantiomorph = 0;

  for (i = 0; i < EndOfAxis; i++) N_count[i] = 0;

  lsmx  = SgInfo->ListSeitzMx;
  lrmxi = SgInfo->ListRotMxInfo;

  if (lrmxi == NULL) lrmxi = &RotMxInfo;

  for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
  {
    if (lrmxi == &RotMxInfo)
    {
      if (GetRotMxInfo(lsmx->s.R, &RotMxInfo) == 0) {
        SetSgError(Err_Ill_SMx_in_List);
        return XS_Unknown;
      }
    }

    switch(abs(lrmxi->Order))
    {
      case 1:  i = N1; break;
      case 2:  i = N2; break;
      case 3:  i = N3; break;
      case 4:  i = N4; break;
      case 6:  i = N6; break;
      default:
        SetSgError("Internal Error: FindXtalSystem(): Corrupt ListRotMxInfo");
        return XS_Unknown;
    }

    if (lrmxi->Inverse == 0) /* skip N^-1 */
      N_count[i]++;

    if (lrmxi != &RotMxInfo) lrmxi++;
  }

  i = EndOfAxis;

  if (SgInfo->InversionOffOrigin == 1)
  {
    for (i = 0; i < EndOfAxis; i++)
    {
      if (N_count[i] % 2) break;
      N_count[i] /= 2;
    }
  }

  if (i == EndOfAxis)
  {
    if      (N_count[N3] == 4) SgInfo->XtalSystem = XS_Cubic;
    else if (N_count[N3] >  1) SgInfo->XtalSystem = XS_Unknown;
    else if (N_count[N6] == 1) SgInfo->XtalSystem = XS_Hexagonal;
    else if (N_count[N6] >  0) SgInfo->XtalSystem = XS_Unknown;
    else if (N_count[N3] == 1) SgInfo->XtalSystem = XS_Trigonal;
    else if (N_count[N4] == 1) SgInfo->XtalSystem = XS_Tetragonal;
    else if (N_count[N4] >  0) SgInfo->XtalSystem = XS_Unknown;
    else if (N_count[N2] >  2) SgInfo->XtalSystem = XS_Orthorhombic;
    else if (N_count[N2] >  0) SgInfo->XtalSystem = XS_Monoclinic;
    else if (N_count[N1] >  0) SgInfo->XtalSystem = XS_Triclinic;

    HonorSign = 1;
    iList = FindSeitzMx(SgInfo, -1, HonorSign, 'o', '.');
    if (iList < 0) HonorSign = 0;

    switch(SgInfo->XtalSystem)
    {
      case XS_Monoclinic:
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '=');
        if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
        break;
      case XS_Tetragonal:
        CheckEnantiomorph = 1;
        iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
        if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
        break;
      case XS_Trigonal:
        CheckEnantiomorph = 1;
          iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '=');
        if (iList < 0)
          iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '*');
        if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
        break;
      case XS_Hexagonal:
        CheckEnantiomorph = 1;
        iList = FindSeitzMx(SgInfo, 6, HonorSign, 0, '=');
        if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
        break;
      case XS_Cubic:
        iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
        if (iList >= 0) CheckEnantiomorph = 1;
        break;
      default:
        iList = -1;
        break;
    }
  }

  if (SgInfo->XtalSystem == XS_Unknown) {
    SetSgError("Error: Can't determine crystal system");
  }
  else if (iList >= 0)
  {
        lrmxi = ListOrBufRotMxInfo(SgInfo, iList, &RotMxInfo);
    if (lrmxi == NULL) {
      SgInfo->XtalSystem = XS_Unknown;
      return XS_Unknown;
    }

    if (SgInfo->XtalSystem != XS_Cubic)
    {
      SgInfo->UniqueRefAxis = lrmxi->RefAxis;
      SgInfo->UniqueDirCode = lrmxi->DirCode;

      if (SgInfo->XtalSystem == XS_Trigonal && lrmxi->DirCode == '=')
      {
        switch (lrmxi->RefAxis)
        {
          case 'z':
            switch (SgInfo->LatticeInfo->Code)
            {
              case 'R': SgInfo->ExtraInfo = EI_Obverse; break;
              case 'T': SgInfo->ExtraInfo = EI_Reverse; break;
            }
            break;
          case 'y':
            switch (SgInfo->LatticeInfo->Code)
            {
              case 'S': SgInfo->ExtraInfo = EI_Obverse; break;
              case 'R': SgInfo->ExtraInfo = EI_Reverse; break;
            }
            break;
          case 'x':
            switch (SgInfo->LatticeInfo->Code)
            {
              case 'T': SgInfo->ExtraInfo = EI_Obverse; break;
              case 'S': SgInfo->ExtraInfo = EI_Reverse; break;
            }
            break;
        }
      }
    }

    if (   HonorSign == 0 /* no inversion matrix */
        && SgInfo->LatticeInfo->Code == 'P'
        && CheckEnantiomorph == 1)
    {
      lsmx = &SgInfo->ListSeitzMx[iList];

      if (GetRotMxOrder(lsmx->s.R) > 1)
      {
        i =  lsmx->s.T[0] * lrmxi->EigenVector[0];
        i += lsmx->s.T[1] * lrmxi->EigenVector[1];
        i += lsmx->s.T[2] * lrmxi->EigenVector[2];

        if (i % (STBF / 2)) SgInfo->ExtraInfo = EI_Enantiomorphic;
      }
    }
  }

  return SgInfo->XtalSystem;
}


static int BuildGenerator_iList(T_SgInfo *SgInfo)
{
  int  iList, iList_1, nG;
  int  SgInfo_CI, HonorSign, Flag3asterisk, FlagPG;


  SgInfo_CI = (SgInfo->Centric || SgInfo->InversionOffOrigin);

  SgInfo->PointGroup = PG_Unknown;

        nG   = SgInfo->nGenerator = 0;
#define G_iL   SgInfo->Generator_iList

  HonorSign = 1;
  iList_1 = FindSeitzMx(SgInfo, -1, HonorSign, 'o', '.');
  if (iList_1 < 0) HonorSign = 0;

  if (SgInfo->XtalSystem == XS_Unknown)
    FindXtalSystem(SgInfo);

  switch(SgInfo->XtalSystem)
  {
    case XS_Triclinic:
      if (iList_1 < 0)
        iList_1 = FindSeitzMx(SgInfo, 1, HonorSign, 'o', '.');
      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      if (SgInfo_CI)
        SgInfo->PointGroup = PG_1b;
      else
        SgInfo->PointGroup = PG_1;

      SgInfo->nGenerator = nG;
      return 0;

    case XS_Monoclinic:
      iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '=');
      if (iList < 0) break;
      G_iL[nG++] = iList;

      if (SgInfo_CI)
        SgInfo->PointGroup = PG_2_m;
      else if (deterRotMx(SgInfo->ListSeitzMx[iList].s.R) == -1)
        SgInfo->PointGroup = PG_m;
      else
        SgInfo->PointGroup = PG_2;

      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      SgInfo->nGenerator = nG;
      return 0;

    case XS_Orthorhombic:
      iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
      if (iList >= 0) G_iL[nG++] = iList;

      iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
      if (iList >= 0) G_iL[nG++] = iList;

      if (nG < 2)
      {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 'y', '=');
        if (iList >= 0) G_iL[nG++] = iList;
      }

      if (nG != 2) break;

      if (SgInfo_CI)
        SgInfo->PointGroup = PG_mmm;
      else if (   deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1
               || deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
        SgInfo->PointGroup = PG_mm2;
      else
        SgInfo->PointGroup = PG_222;

      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      SgInfo->nGenerator = nG;
      return 0;

    case XS_Tetragonal:
      iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
      if (iList < 0) break;
      G_iL[nG++] = iList;

      if (          SgInfo->UniqueRefAxis != 'x')
      {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
        if (iList >= 0) G_iL[nG++] = iList;
      }
      if (nG < 2 && SgInfo->UniqueRefAxis != 'z')
      {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
        if (iList >= 0) G_iL[nG++] = iList;
      }
      if (nG < 2 && SgInfo->UniqueRefAxis != 'y')
      {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 'y', '=');
        if (iList >= 0) G_iL[nG++] = iList;
      }

      if (nG < 2)
      {
        if (SgInfo_CI)
          SgInfo->PointGroup = PG_4_m;
        else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
          SgInfo->PointGroup = PG_4b;
        else
          SgInfo->PointGroup = PG_4;
      }
      else
      {
        if (SgInfo_CI)
          SgInfo->PointGroup = PG_4_mmm;
        else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
        {
          if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
            SgInfo->PointGroup = PG_4bm2;
          else
            SgInfo->PointGroup = PG_4b2m;
        }
        else
        {
          if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
            SgInfo->PointGroup = PG_4mm;
          else
            SgInfo->PointGroup = PG_422;
        }
      }

      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      SgInfo->nGenerator = nG;
      return 0;

    case XS_Trigonal:
    case XS_Hexagonal:
      Flag3asterisk = 0;

      if (SgInfo->XtalSystem == XS_Trigonal)
      {
          iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '=');
        if (iList < 0)
        {
          iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '*');
          Flag3asterisk = 1;
        }
      }
      else
        iList = FindSeitzMx(SgInfo, 6, HonorSign, 0, '=');

      if (iList < 0) break;
      G_iL[nG++] = iList;

      iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '\'');
      if (iList >= 0) G_iL[nG++] = iList;

      FlagPG = -1;

      if (nG < 2)
      {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '"');
        if (iList >= 0) G_iL[nG++] = iList;
        FlagPG = 1;
      }

      if (SgInfo->XtalSystem == XS_Trigonal)
      {
        if (nG < 2)
        {
          if (SgInfo_CI) SgInfo->PointGroup = PG_3b;
          else           SgInfo->PointGroup = PG_3;
        }
        else
        {
          if (Flag3asterisk == 1)
          {
            if (SgInfo_CI)
              SgInfo->PointGroup = PG_3bm;
            else
            {     FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
              if (FlagPG == -1) SgInfo->PointGroup = PG_3m;
              else              SgInfo->PointGroup = PG_32;
            }
          }
          else if (FlagPG == -1)
          {
            if (SgInfo_CI)
              SgInfo->PointGroup = PG_3b1m;
            else
            {     FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
              if (FlagPG == -1) SgInfo->PointGroup = PG_31m;
              else              SgInfo->PointGroup = PG_312;
            }
          }
          else
          {
            if (SgInfo_CI)
              SgInfo->PointGroup = PG_3bm1;
            else
            {     FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
              if (FlagPG == -1) SgInfo->PointGroup = PG_3m1;
              else              SgInfo->PointGroup = PG_321;
            }
          }
        }
      }
      else
      {
        if (nG < 2)
        {
          if (SgInfo_CI)
            SgInfo->PointGroup = PG_6_m;
          else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
            SgInfo->PointGroup = PG_6b;
          else
            SgInfo->PointGroup = PG_6;
        }
        else
        {
          if (SgInfo_CI)
            SgInfo->PointGroup = PG_6_mmm;
          else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
          {
            if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == FlagPG)
              SgInfo->PointGroup = PG_6b2m;
            else
              SgInfo->PointGroup = PG_6bm2;
          }
          else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
            SgInfo->PointGroup = PG_6mm;
          else
            SgInfo->PointGroup = PG_622;
        }
      }

      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      SgInfo->nGenerator = nG;
      return 0;

    case XS_Cubic:
      FlagPG = 0;

        iList = FindSeitzMx(SgInfo, 4, HonorSign, 'z', '=');
      if (iList < 0) {
        iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
        FlagPG = 1;
      }
      if (iList < 0) break;
      G_iL[nG++] = iList;

      iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
      if (iList < 0) break;
      G_iL[nG++] = iList;

      iList = FindSeitzMx(SgInfo, 3, HonorSign, 'o', '*');
      if (iList < 0) break;
      G_iL[nG++] = iList;

      if (FlagPG)
      {
        if (SgInfo_CI) SgInfo->PointGroup = PG_m3b;
        else           SgInfo->PointGroup = PG_23;
      }
      else
      {
        if (SgInfo_CI)
          SgInfo->PointGroup = PG_m3bm;
        else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
          SgInfo->PointGroup = PG_4b3m;
        else
          SgInfo->PointGroup = PG_432;
      }

      if (iList_1 >= 0) G_iL[nG++] = iList_1;

      SgInfo->nGenerator = nG;
      return 0;

    default:
      break;
  }

#undef G_iL

  return -1;
}


static int BuildHSym(T_SgInfo *SgInfo)
{
  int                NeedDash, HaveOrSh, nGRT, iList, iG, ip, os, i;
  int                AbsOrder, RefAxis, DirCode, ScrewPrimitive, Screw;
  int                PreviousRotation, PreviousRefAxis, PreviousDirCode;
  int                nTrV, iTrV, OrSh[3], RO[3], Transl[3];
  const int          *TrV, *ht;
  T_RTMx             SMx_1;
  const T_RTMx       *SMx;
  const T_RotMxInfo  *rmxi;
  char               *hsym, *hsym_mark;

  struct
  {
    T_RotMxInfo         RMxI_Buf;
    const T_RotMxInfo  *RMxI;
    int                Transl[3];
  }
  GRT[sizeof SgInfo->Generator_iList / sizeof (*SgInfo->Generator_iList) + 1];

  const char *Digits = "0123456";


  if (SgInfo->nGenerator == 0) {
    SetSgError("Internal Error: BuildHSym(): Empty generator list");
    return -1;
  }

  HaveOrSh = 0;

  for (i = 0; i < 3; i++) {
        OrSh[i] = SgInfo->OriginShift[i] * (STBF / 12);
    if (OrSh[i]) HaveOrSh = 1;
  }

  NeedDash = 0;
  nGRT = 0;

  for (iG = 0; iG < SgInfo->nGenerator; iG++)
  {
    iList = SgInfo->Generator_iList[iG];

        GRT[nGRT].RMxI = ListOrBufRotMxInfo(SgInfo, iList, &GRT[nGRT].RMxI_Buf);
    if (GRT[nGRT].RMxI == NULL)
      return -1;

    SMx = &SgInfo->ListSeitzMx[iList];

    RotMx_t_Vector(RO, SMx->s.R, OrSh, STBF);

    for (i = 0; i < 3; i++)
      GRT[nGRT].Transl[i] = iModPositive(SMx->s.T[i] + RO[i] - OrSh[i], STBF);

    if (GRT[nGRT].RMxI->Order == -1)
    {
      for (i = 0; i < 3; i++)
        if (GRT[nGRT].Transl[i] != 0) break;

      if (i == 3) NeedDash = 1;
      else        nGRT++;
    }
    else
      nGRT++;
  }

  if (SgInfo->Centric)
  {
    if (HaveOrSh == 0)
      NeedDash = 1;
    else
    {
      for (iG = 0; iG < nGRT; iG++)
        if (GRT[iG].RMxI->Order == 1) break;

      InitSeitzMx(&SMx_1, -1);

      if (GetRotMxInfo(SMx_1.s.R, &GRT[iG].RMxI_Buf) != -1) {
        SetSgError("Internal Error: BuildHSym(): Corrupt GetRotMxInfo()");
        return -1;
      }

      GRT[iG].RMxI = &GRT[iG].RMxI_Buf;

      for (i = 0; i < 3; i++)
        GRT[iG].Transl[i] = iModPositive(-2 * OrSh[i], STBF);

      if (iG == nGRT)
        nGRT++;
    }
  }

  hsym = SgInfo->HallSymbol;

  for (i = 0; i <= MaxLenHallSymbol; i++)
    *hsym++ = '\0';

  PreviousRotation = 0;
  PreviousRefAxis = 0;
  PreviousDirCode = 0;

  hsym = SgInfo->HallSymbol;

  if (NeedDash)
    *hsym++ = '-';
  else
    *hsym++ = ' ';

  *hsym++ = SgInfo->LatticeInfo->Code;

  nTrV = SgInfo->LatticeInfo->nTrVector;

  for (iG = 0; iG < nGRT; iG++)
  {
    rmxi = GRT[iG].RMxI;

    AbsOrder = abs(rmxi->Order);
    RefAxis = rmxi->RefAxis;
    DirCode = rmxi->DirCode;
    if (RefAxis == 'o') RefAxis = 0;
    if (DirCode == '=' || DirCode == '.') DirCode = 0;

    if (iG == 0)
    {
      if (RefAxis == 'z') RefAxis = 0;
    }
    else
    {
      if      (AbsOrder == 2)
      {
        if      (PreviousRotation == 2 || PreviousRotation == 4)
        {
          if (RefAxis == 'x') RefAxis = 0;
        }
        else if (PreviousRotation == 3 || PreviousRotation == 6)
        {
          if (   PreviousDirCode == '*'
              || RefAxis == PreviousRefAxis) RefAxis = 0;
          if (DirCode == '\'') DirCode = 0;
        }
      }
      else if (AbsOrder == 3)
      {
        if (DirCode == '*') DirCode = 0;
      }
    }

    PreviousRotation = AbsOrder;
    PreviousRefAxis = rmxi->RefAxis;
    PreviousDirCode = rmxi->DirCode;

    *hsym++ = ' ';
    if (rmxi->Order < 0) *hsym++ = '-';
    *hsym++ = Digits[AbsOrder];
    if (RefAxis) *hsym++ = RefAxis;
    if (DirCode) *hsym++ = DirCode;

    TrV = SgInfo->LatticeInfo->TrVector;

    for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
    {
      for (i = 0; i < 3; i++)
        if ((GRT[iG].Transl[i] + TrV[i]) % STBF != 0)
          break;

      if (i == 3)
        break;
    }

    if (iTrV < nTrV)
      continue; /* next iG */

    hsym_mark = hsym;

    TrV = SgInfo->LatticeInfo->TrVector;

    for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3, hsym = hsym_mark)
    {
      for (i = 0; i < 3; i++)
        Transl[i] = iModPositive(GRT[iG].Transl[i] + TrV[i], STBF);

      Screw = 0;

      for (ip = 0; ip < 3; ip++)
        if (rmxi->EigenVector[ip] != 0) break;

      if (ip < 3 && rmxi->EigenVector[ip] == 1)
      {
        for (i = ip + 1; i < 3; i++)
          if (rmxi->EigenVector[i] != 0) break;

        if (i == 3)
        {
          ScrewPrimitive = STBF / AbsOrder;
          Screw = Transl[ip] / ScrewPrimitive;
              i = Screw * ScrewPrimitive;
          if (i % 3)
          {
            *hsym++ = Digits[Screw];
            Transl[ip] -= i;
          }
        }
      }

      ht = HallTranslations;

      while (*ht)
      {
        for (i = 0; i < 3; i++)
          if (Transl[i] < ht[i + 1]) break;

        if (i == 3)
        {
          for (i = 0; i < 3; i++)
            Transl[i] -= ht[i + 1];

          *hsym++ = (char) *ht;
        }

        ht += 4;
      }

      for (i = 0; i < 3; i++)
        if (Transl[i] != 0)
          break;

      if (i == 3)
        break;
    }

    if (iTrV == nTrV)
      return 0;
  }

  if (nGRT == 0)
  {
    *hsym++ = ' ';
    *hsym++ = '1';
  }

  if (HaveOrSh)
  {
    *hsym++ = ' ';
    *hsym++ = '(';

    for (i = 0; i < 3; i++)
    {
      if (i) *hsym++ = ' ';

          os = iModPositive(SgInfo->OriginShift[i], 12);
      if (os > 6)
      {
        *hsym++ = '-';
        os = 12 - os;
      }

      *hsym++ = Digits[os];
    }

    *hsym++ = ')';
  }

  *hsym = '\0';

  if (SgInfo->HallSymbol[MaxLenHallSymbol] != '\0') {
    SetSgError("Internal Error: BuildHSym(): MaxLenHallSymbol too small");
    return -1;
  }

  return 1;
}


static int BuildHallSymbol(T_SgInfo *SgInfo, int FixedOriginShift)
{
  int     ix, iy, iz;
  int     status;

  static const int ShiftTable[] = { 0, 1, -1, 2, -2, 3 };


  if (SgError != NULL) return -1;

  if (SgInfo->nGenerator == 0)
  {
    if (BuildGenerator_iList(SgInfo) != 0)
    {
      SetSgError("Error: Can't build generator list");
      return -1;
    }
  }

  if (FixedOriginShift)
  {
        status = BuildHSym(SgInfo);
    if (status == 1)
      return 0;
  }
  else
  {
    for (ix = 0; ix < 6; ix++)
    {
      SgInfo->OriginShift[0] = ShiftTable[ix];

      for (iy = 0; iy < 6; iy++)
      {
        SgInfo->OriginShift[1] = ShiftTable[iy];

        for (iz = 0; iz < 6; iz++)
        {
          SgInfo->OriginShift[2] = ShiftTable[iz];

              status = BuildHSym(SgInfo);
          if (status < 0)
            return -1;

          if (status == 1)
            return 0;
        }
      }
    }
  }

  SetSgError("Error: Can't build Hall Symbol");
  return -1;
}


void InitSgInfo(T_SgInfo *SgInfo)
{
  int  i;


  SgInfo->GenOption = 0;
  SgInfo->Centric = 0;
  SgInfo->InversionOffOrigin = 0;
  SgInfo->LatticeInfo = LI_P;
  SgInfo->StatusLatticeTr = 0;
  for (i = 0; i < 3; i++)
    SgInfo->OriginShift[i] = 0;
  SgInfo->nList = 0;

  SgInfo->OrderL = 0;
  SgInfo->OrderP = 0;
  SgInfo->XtalSystem = XS_Unknown;
  SgInfo->UniqueRefAxis = 0;
  SgInfo->UniqueDirCode = 0;
  SgInfo->ExtraInfo = EI_Unknown;
  SgInfo->PointGroup = PG_Unknown;
  SgInfo->nGenerator = 0;
  SgInfo->HallSymbol[0] = '\0';
  SgInfo->TabSgName = NULL;
  SgInfo->CCMx_LP = NULL;
  SgInfo->n_si_Vector = -1;
}


int CompleteSgInfo(T_SgInfo *SgInfo)
{
  int                List_iList[192];
  const T_TabSgName  *tsgn;


  if (SgInfo->StatusLatticeTr == -1)
  {
    if (AddLatticeTr2ListSeitzMx(SgInfo, SgInfo->LatticeInfo) < 0)
      return -1;
  }

  if (ApplyOriginShift(SgInfo) < 0)
    return -1;

  if (SgInfo->nList > sizeof List_iList / sizeof (*List_iList)) {
    SetSgError("Internal Error: CompleteSgInfo()");
    return -1;
  }

  if (SgInfo->nList > 1)
  {
    SortSgInfoList(SgInfo, List_iList);
    if (SgError != NULL) return -1;
  }

  if (RemoveLatticeTr(SgInfo) != 0)
    return -1;

  if (RemoveInversion(SgInfo) != 0)
    return -1;

  TidyTranslation(SgInfo);

  if (SgInfo->nList > 1)
  {
    SortSgInfoList(SgInfo, List_iList);
    if (SgError != NULL) return -1;
  }
                             SgInfo->OrderP = SgInfo->nList;
  if (SgInfo->Centric == -1) SgInfo->OrderP *= 2;

  SgInfo->OrderL = SgInfo->OrderP * SgInfo->LatticeInfo->nTrVector;

  if (BuildHallSymbol(SgInfo, 0) != 0)
    return -1;

  for (tsgn = TabSgName; tsgn->HallSymbol; tsgn++)
    if (   strcmp(tsgn->HallSymbol, SgInfo->HallSymbol) == 0
        && (   SgInfo->TabSgName == NULL
            || SgInfo->TabSgName == tsgn))
      break;

  if (SgInfo->TabSgName != NULL && tsgn->HallSymbol == NULL)
  {
    if (SgError) return -1;

    sprintf(SgErrorBuffer,
      "Internal Error: Input/Output HallSymbol mismatch: %s <> %s",
      SgInfo->TabSgName->HallSymbol, SgInfo->HallSymbol);

    SetSgError(SgErrorBuffer);
    return -1;
  }

  if (tsgn->HallSymbol)
    SgInfo->TabSgName = tsgn;

  SgInfo->CCMx_LP = NULL;

  switch (SgInfo->LatticeInfo->Code)
  {
    case 'P': SgInfo->CCMx_LP = CCMx_PP; break;
    case 'A': SgInfo->CCMx_LP = CCMx_AP; break;
    case 'B': SgInfo->CCMx_LP = CCMx_BP; break;
    case 'C': SgInfo->CCMx_LP = CCMx_CP; break;
    case 'I': SgInfo->CCMx_LP = CCMx_IP; break;
    case 'R':
      switch (SgInfo->UniqueRefAxis) {
        case   0:
        case 'z': SgInfo->CCMx_LP = CCMx_RP_z; break;
        case 'y': SgInfo->CCMx_LP = CCMx_RP_y; break;
        default: break;
      }
      break;
    case 'S':
      switch (SgInfo->UniqueRefAxis) {
        case   0:
        case 'y': SgInfo->CCMx_LP = CCMx_SP_y; break;
        case 'x': SgInfo->CCMx_LP = CCMx_SP_x; break;
        default: break;
      }
      break;
    case 'T':
      switch (SgInfo->UniqueRefAxis) {
        case   0:
        case 'x': SgInfo->CCMx_LP = CCMx_TP_x; break;
        case 'z': SgInfo->CCMx_LP = CCMx_TP_z; break;
        default: break;
      }
      break;
    case 'F': SgInfo->CCMx_LP = CCMx_FP; break;
    default:
      break;
  }

  if (SgInfo->CCMx_LP == NULL) {
    SetSgError("Internal Error: Illegal lattice code");
    return -1;
  }

  return 0;
}


int CB_SMx(T_RTMx *CSiC,
           const T_RTMx *CBMx, const T_RTMx *SMx, const T_RTMx *InvCBMx)
{
  int     i;
  T_RTMx  BufMx;


  RTMxMultiply(&BufMx, SMx,  InvCBMx, CTBF / STBF, CTBF);
  RTMxMultiply(CSiC,   CBMx, &BufMx,  CRBF,        CRBF * CTBF);

  for (i = 0; i < 9; i++)
  {
    if (CSiC->s.R[i] % (CRBF * CRBF)) {
      SetSgError("Internal Error: Corrupt CBMx/SMx/InvCBMx");
      return -1;
    }

    CSiC->s.R[i] /= (CRBF * CRBF);
  }

  for (i = 0; i < 3; i++)
  {
    if (CSiC->s.T[i] % (CRBF * (CTBF / STBF))) {
      SetSgError("Internal Error: Out of STBF range");
      return -1;
    }

    CSiC->s.T[i] /= (CRBF * (CTBF / STBF));
  }

  return 0;
}


int TransformSgInfo(const T_SgInfo *SgInfo,
                    const T_RTMx *CBMx, const T_RTMx *InvCBMx,
                    T_SgInfo *BC_SgInfo)
{
  int           iList, f, i;
  int           nTrV, iTrV, nLoopInv, iLoopInv;
  const int     *TrV;
  T_RTMx        SMx, BC_SMx;
  const T_RTMx  *lsmx;


  nLoopInv = Sg_nLoopInv(SgInfo);

  nTrV = SgInfo->LatticeInfo->nTrVector;
   TrV = SgInfo->LatticeInfo->TrVector;

  for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
  {
    for (iLoopInv = 0; iLoopInv < nLoopInv; iLoopInv++)
    {
      if (iLoopInv == 0) f =  1;
      else               f = -1;

      lsmx = SgInfo->ListSeitzMx;

      for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
      {
        for (i = 0; i < 9; i++)
          SMx.s.R[i] = f * lsmx->s.R[i];

        for (i = 0; i < 3; i++)
          SMx.s.T[i] = f * lsmx->s.T[i] + TrV[i];

        if (CB_SMx(&BC_SMx, CBMx, &SMx, InvCBMx) != 0)
          return -1;

        if (Add2ListSeitzMx(BC_SgInfo, &BC_SMx) < 0)
          return -1;
      }
    }
  }

  return 0;
}


#undef SGCLIB_C__
