sgmath.c

/* $Id: sgmath.c 851 2002-02-09 05:39:17Z wdelano $ */

/* The source code contained in this file is            */
/* Copyright (C) 1994-2000 by Ralf W. Grosse-Kunstleve. */
/* Please see the LICENSE file for more information.    */

#ifdef _PYMOL_WIN32
#include"os_predef.h"
#endif

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


#undef SG_GLOBAL
#include "sglite.h"


int iGCD(const int a, const int b)
{
  int  ri, rj, rk;


      ri = a;
  if (ri < 0) ri = -ri;

  if ((rj = b) != 0)
  {
    for (;;)
    {
      rk = ri % rj; if (rk == 0) { ri = rj; break; }
      ri = rj % rk; if (ri == 0) { ri = rk; break; }
      rj = rk % ri; if (rj == 0) {          break; }
    }

    if (ri < 0) ri = -ri;
  }

  return ri;
}


int FindGCD(const int *S, int nS)
{
  int  ri, rj, rk;


  if (nS-- == 0) return 0;

      ri = *S++;
  if (ri < 0) ri = -ri;

  while (nS--)
  {
    if ((rj = *S++) != 0)
    {
      for (;;)
      {
        rk = ri % rj; if (rk == 0) { ri = rj; break; }
        ri = rj % rk; if (ri == 0) { ri = rk; break; }
        rj = rk % ri; if (rj == 0) {          break; }
      }

      if (ri < 0) ri = -ri;

      if (ri == 1) break;
    }
  }

  return ri;
}


int CancelGCD(int *S, int nS)
{
  int  GCD, i;


      GCD = FindGCD(S, nS);
  if (GCD)
    rangei(3) S[i] /= GCD;

  return GCD;
}


int CancelBFGCD(int *S, int nS, int BF)
{
  int  GCD, i;

      GCD = FindGCD(S, nS);
      GCD = iGCD(GCD, BF);
  if (GCD) {
    rangei(3) S[i] /= GCD;
    return BF / GCD;
  }

  return GCD;
}


int iLCM(const int a, const int b)
{
  int  al, ri, rj, rk;


      al = a;
  if (al == 0) al = 1;

       ri = al;
  if ((rj = b) != 0)
  {
    for (;;)
    {
      rk = ri % rj; if (rk == 0) { ri = rj; break; }
      ri = rj % rk; if (ri == 0) { ri = rk; break; }
      rj = rk % ri; if (rj == 0) {          break; }
    }

    ri = al / ri * b;
  }

  if (ri < 0) return -ri;
              return  ri;
}


int FindLCM(const int *S, int nS)
{
  int  a, b, ri, rj, rk;


  if (nS-- == 0) return 1;

      ri = *S++;
  if (ri == 0) ri = 1;

  while (nS--)
  {
    if ((rj = *S++) != 0)
    {
      a = ri; b = rj;

      for (;;)
      {
        rk = ri % rj; if (rk == 0) { ri = rj; break; }
        ri = rj % rk; if (ri == 0) { ri = rk; break; }
        rj = rk % ri; if (rj == 0) {          break; }
      }

      ri = a / ri * b;
    }
  }

  if (ri < 0) return -ri;
              return  ri;
}


void SimplifyFraction(int nume, int deno, int *o_nume, int *o_deno)
{
  int gcd = iGCD(nume, deno);
  if (gcd)
  {
    *o_nume = nume / gcd;
    *o_deno = deno / gcd;

    if (*o_deno < 0) {
      *o_nume *= -1;
      *o_deno *= -1;
    }
  }
}


#define GmulV(GV, G, V) \
  (GV)[0] = (G)[0] * (V)[0] + (G)[1] * (V)[1] + (G)[2] * (V)[2]; \
  (GV)[1] = (G)[3] * (V)[0] + (G)[4] * (V)[1] + (G)[5] * (V)[2]; \
  (GV)[2] = (G)[6] * (V)[0] + (G)[7] * (V)[1] + (G)[8] * (V)[2]


int iScalProd(const int *u, const int *v, const int *PseudoG)
{
  int  Prod, Gv[3];


  if (PseudoG)
  {
    GmulV(Gv, PseudoG, v);
    v = Gv;
  }

  Prod =   u[0] * v[0]
         + u[1] * v[1]
         + u[2] * v[2];

  return Prod;
}


void iCrossProd(int *rxs, const int *r, const int *s, const int *PseudoG)
{
  int  Gr[3], Gs[3];


  if (PseudoG)
  {
    GmulV(Gr, PseudoG, r);
    GmulV(Gs, PseudoG, s);
    r = Gr;
    s = Gs;
  }

  rxs[0] = r[1] * s[2] - r[2] * s[1];
  rxs[1] = r[2] * s[0] - r[0] * s[2];
  rxs[2] = r[0] * s[1] - r[1] * s[0];
}


#undef GmulV


int AreLinDepV(const int a[3], const int b[3])
{
  int        axb[3], i;
  const int  V000[3] = { 0, 0, 0 };


  iCrossProd(axb, a, b, NULL);
  if (MemCmp(axb, V000, 3) == 0)
  {
    rangei(3) {
      if (a[i]) {
        if (abs(a[i]) > abs(b[i])) return 1;
        return -1;
      }
    }
  }

  return 0;
}


void SMx_t_InvT(const T_RTMx *SMx, const int InvT[3], T_RTMx *ProdSMx)
{
  int   i;

  rangei(9) ProdSMx->s.R[i] = -SMx->s.R[i];
  rangei(3) ProdSMx->s.T[i] = -SMx->s.T[i] + InvT[i];
}


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;
  }
}


int CBMxMultiply(T_RTMx *ab, const T_RTMx *a, const T_RTMx *b)
{
  T_RTMx  abf[1];

  RTMxMultiply(abf, a, b, CRBF, CRBF * CTBF);
  if (ChangeBaseFactor(abf->a, CRBF, ab->a, 1, 12) != 0)
    return IE(-1);

  return 0;
}


int CBMx2Multiply(T_RTMx ab[2], const T_RTMx a[2], const T_RTMx b[2])
{
  if (CBMxMultiply(&ab[0], &a[0], &b[0]) != 0) return -1;
  if (CBMxMultiply(&ab[1], &b[1], &a[1]) != 0) return -1;
  return 0;
}


int CBMx2Update(T_RTMx ab[2], const T_RTMx a[2])
{
  if (CBMxMultiply(&ab[0],  &a[0], &ab[0]) != 0) return -1;
  if (CBMxMultiply(&ab[1], &ab[1],  &a[1]) != 0) return -1;
  return 0;
}


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 iCoFactorMxTp(const int *Mx, int *CFMxTp)
{
  CFMxTp[0] =  Mx[4] * Mx[8] - Mx[5] * Mx[7];
  CFMxTp[1] = -Mx[1] * Mx[8] + Mx[2] * Mx[7];
  CFMxTp[2] =  Mx[1] * Mx[5] - Mx[2] * Mx[4];
  CFMxTp[3] = -Mx[3] * Mx[8] + Mx[5] * Mx[6];
  CFMxTp[4] =  Mx[0] * Mx[8] - Mx[2] * Mx[6];
  CFMxTp[5] = -Mx[0] * Mx[5] + Mx[2] * Mx[3];
  CFMxTp[6] =  Mx[3] * Mx[7] - Mx[4] * Mx[6];
  CFMxTp[7] = -Mx[0] * Mx[7] + Mx[1] * Mx[6];
  CFMxTp[8] =  Mx[0] * Mx[4] - Mx[1] * Mx[3];
}


int InverseRotMx(const int R[9], int InvR[9], const int RBF)
{
  /*
     InvR = Inverse[R]
     DetF = |R| * RBF^3
   */

  int  DetF, i;

      DetF = deterRotMx(R);
  if (DetF == 0) return 0;

  iCoFactorMxTp(R, InvR);

  rangei(9) InvR[i] *= (RBF * RBF);

  rangei(9) {
    if (InvR[i] %  DetF) return 0;
        InvR[i] /= DetF;
  }

  return DetF;
}


int InverseRTMx(const T_RTMx *RTMx, T_RTMx *InvRTMx, int RBF)
{
  /*
     InvRTMx->s.R =  Inverse[RTMx->s.R]
     InvRTMx->s.T = -Inverse[RTMx->s.R] * RTMx->s.T
     DetF = |RTMx| * RBF^3;
   */

  int  DetF, i;

      DetF = InverseRotMx(RTMx->s.R, InvRTMx->s.R, RBF);
  if (DetF == 0) return 0;

  RotMx_t_Vector(InvRTMx->s.T, InvRTMx->s.R, RTMx->s.T, 0);

  for (i = 0; i < 3; i++) {
    if (InvRTMx->s.T[i] %  RBF) return 0;
        InvRTMx->s.T[i] /= RBF;
        InvRTMx->s.T[i] *= -1;
  }

  return DetF;
}


void iMxMultiply(int *ab, const int *a, const int *b,
                 const int ma, const int na, const int nb)
{
  int        i, j, k;
  const int  *ai, *aij, *bk, *bkj;

  ai = a;

  for (i = 0; i < ma; i++)
  {
    bk = b;

    for (k = 0; k < nb; k++)
    {
      aij = ai;
      bkj = bk;

      *ab = 0;

      for (j = 0; j < na; j++)
      {
        *ab += (*aij) * (*bkj);

        aij++;
        bkj += nb;
      }

      ab++;
      bk++;
    }

    ai += na;
  }
}


int *IdentityMat(int *M, int m)
{
  int  i;

  rangei(m * m) M[i] = 0;
  rangei(m)     M[i * (m + 1)] = 1;

  return M;
}


int *TransposedMat(int *M, int mr, int mc)
{
  int  ir, ic, i;
  int  *Mt;


  if (mc <= 0 || mr <= 0) return NULL;

  nxs_malloc(Mt, mc * mr);
  if (Mt == NULL) {
    (void) SetSgNotEnoughCore(0);
    return NULL;
  }
  i = 0; range1(ir, mr) range1(ic, mc) Mt[ic * mr + ir] = M[i++];
  (void) MemCpy(M, Mt, mc * mr);
  free(Mt);

  return M;
}


int iRowEchelonFormT(int *M, int mr, int mc, int *T, int tc)
{
  /* C version of RowEchelonFormT from the CrystGAP package
         (GAP Version 3.4.4).
     B. Eick, F. Ga"hler and W. Nickel
         Computing Maximal Subgroups and Wyckoff Positions of Space Groups
         Acta Cryst. (1997). A53, 467 - 474
   */

  int  a, i, j, k, ic, Cleared;


#define  M2D(i, j) M[i * mc + j]
#define  T2D(i, j) T[i * tc + j]

  for (i = j = 0; i < mr && j < mc;)
  {
    k = i; while (k < mr && M2D(k, j) == 0) k++;

    if (k == mr)
      j++;
    else
    {
      if (i != k) {
               (void) IntSwap(&M2D(i, 0), &M2D(k, 0), mc);
        if (T) (void) IntSwap(&T2D(i, 0), &T2D(k, 0), tc);
      }

      range2(k, k + 1, mr) {
        a = abs(M2D(k, j));
        if (a != 0 && a < abs(M2D(i, j))) {
                 (void) IntSwap(&M2D(i, 0), &M2D(k, 0), mc);
          if (T) (void) IntSwap(&T2D(i, 0), &T2D(k, 0), tc);
        }
      }

      if (M2D(i, j) < 0) {
               range1(ic, mc) M2D(i, ic) *= -1;
        if (T) range1(ic, tc) T2D(i, ic) *= -1;
      }

      Cleared = 1;
      range2(k, i + 1, mr) {
        a = M2D(k, j) / M2D(i, j);
        if (a != 0) {
                 range1(ic, mc) M2D(k, ic) -= a * M2D(i, ic);
          if (T) range1(ic, tc) T2D(k, ic) -= a * T2D(i, ic);
        }
        if (M2D(k, j) != 0) Cleared = 0;
      }
      if (Cleared) { i++; j++; }
    }
  }

#undef  M2D
#undef  T2D

  return i;
}


static int IsDiagonalMat(const int *M, int mr, int mc)
{
  int  ir, ic;

  if (mr != mc) return 0;
  range1(ir, mr) range1(ic, mc) if (ir != ic && M[ir * mc + ic]) return 0;

  return 1;
}


int iREBacksubst(const int *M, const int *V,
                 const int nr, const int nc,
                 int *Sol, int *FlagIndep)
{
  int  ir, ic, icp, nv;
  int  d, m, f, jc;

#define M2D(i, j) M[i * nc + j]

  if (FlagIndep)
    for (ic = 0; ic < nc; ic++) FlagIndep[ic] = 1;

  d = 1;

  for (ir = nr - 1; ir >= 0; ir--)
  {
    range1(ic, nc) if (M2D(ir, ic)) goto Set_Sol_ic;

    if (V && V[ir] != 0) return 0;
    continue;

    Set_Sol_ic:

    if (FlagIndep) FlagIndep[ic] = 0;

    if (Sol) {
                    icp = ic + 1;
          nv = nc - icp;
      if (nv) {
        iMxMultiply(&Sol[ic], &M2D(ir, icp), &Sol[icp], 1, nv, 1);
                     Sol[ic] *= -1;
      }
      else
        Sol[ic] = 0;

      if (V) Sol[ic] += d * V[ir];

                        m = M2D(ir, ic);
      f = iGCD(Sol[ic], m);
      if (m < 0) f *= -1;
      Sol[ic] /= f;
      f = m / f;
      if (f != 1) {
        range1(jc, nc) if (jc != ic) Sol[jc] *= f;
        d *= f;
      }
    }
  }

#undef M2D

  return d;
}


int iRESetIxIndep(const int *REMx, int nr, int nc, int *IxIndep, int mIndep)
{
  int  nIndep, ic;
  int  FlagIndep[6];

  if (nc > sizeof FlagIndep / sizeof (*FlagIndep))
    return IE(-1);

  if (iREBacksubst(REMx, NULL, nr, nc, NULL, FlagIndep) < 1)
    return IE(-1);

  nIndep = 0;

  range1(ic, nc) {
    if (FlagIndep[ic]) {
      if (nIndep == mIndep) return -1;
      IxIndep[nIndep++] = ic;
    }
  }

  return nIndep;
}


int SolveHomRE2(const int REMx[9], int EV[3])
{
  /* REMx must be in row echelon form with Rank 2.
   */

  int  IxIndep[1], i;

  if (iRESetIxIndep(REMx, 2, 3, IxIndep, 1) != 1)
    return IE(-1);

  rangei(3) EV[i] = 0;
  EV[IxIndep[0]] = 1;

  if (iREBacksubst(REMx, NULL, 2, 3, EV, NULL) < 1)
    return IE(-1);

  if (SignHemisphere(EV[0], EV[1], EV[2]) < 0)
    rangei(3) EV[i] *= -1;

  return 0;
}


int SolveHomRE1(const int REMx[3], const int IxIndep[2], int Sol[4][3])
{
  int        iPV, i;
  const int  TrialV[4][2] =
    {{ 1,  0 },
     { 0,  1 },
     { 1,  1 },
     { 1, -1 },
    };

  range1(iPV, 4)
  {
    rangei(3) Sol[iPV][i] = 0;
    rangei(2) Sol[iPV][IxIndep[i]] = TrialV[iPV][i];

    if (iREBacksubst(REMx, NULL, 2, 3, Sol[iPV], NULL) < 1)
      return IE(-1);
  }

  return 0;
}


int SmithNormalForm(int *M, int mr, int mc, int *P, int *Q)
{
  int  rr, rc;


  rr = mr;
  rc = mc;

  if (P) (void) IdentityMat(P, mr);
  if (Q) (void) IdentityMat(Q, mc);

  for (;;)
  {
    rr = iRowEchelonFormT(M, rr, rc, P, mr);
    if (IsDiagonalMat(M, rr, rc)) break;
    (void) TransposedMat(M, rr, rc);

    rc = iRowEchelonFormT(M, rc, rr, Q, mc);
    if (IsDiagonalMat(M, rc, rr)) break;
    (void) TransposedMat(M, rc, rr);
  }

  return rr;
}