/*
  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__