text: /* Copyright (C) 2010 Atsushi Togo */
  /* All rights reserved. */
  
  /* This file is part of spglib. */
  
  /* Redistribution and use in source and binary forms, with or without */
  /* modification, are permitted provided that the following conditions */
  /* are met: */
  
  /* * Redistributions of source code must retain the above copyright */
  /*   notice, this list of conditions and the following disclaimer. */
  
  /* * Redistributions in binary form must reproduce the above copyright */
  /*   notice, this list of conditions and the following disclaimer in */
  /*   the documentation and/or other materials provided with the */
  /*   distribution. */
  
  /* * Neither the name of the spglib project nor the names of its */
  /*   contributors may be used to endorse or promote products derived */
  /*   from this software without specific prior written permission. */
  
  /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
  /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
  /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
  /* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
  /* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
  /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
  /* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
  /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
  /* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
  /* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
  /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
  /* POSSIBILITY OF SUCH DAMAGE. */
  
  #include <stdio.h>
  #include "debug.h"
  #include "hall_symbol.h"
  #include "spg_database.h"
  #include "spacegroup.h"
  #include "symmetry.h"
  #include "mathfunc.h"
  
  static int M_bcc[3][3] = {
    { 0, 1, 1 },
    { 1, 0, 1 },
    { 1, 1, 0 },
  };
  static int M_fcc[3][3] = {
    {-1, 1, 1 },
    { 1,-1, 1 },
    { 1, 1,-1 },
  };
  static int M_ac[3][3] = {
    { 1, 0, 0 },
    { 0, 1, 1 },
    { 0,-1, 1 },
  };
  static int M_bc[3][3] = {
    { 1, 0, 1 },
    { 0, 1, 0 },
    {-1, 0, 1 },
  };
  static int M_cc[3][3] = {
    { 1,-1, 0 },
    { 1, 1, 0 },
    { 0, 0, 1 },
  };
  static int M_rc[3][3] = {
    { 1, 0, 1 },
    {-1, 1, 1 },
    { 0,-1, 1 },
  };
  
  static double M_bcc_inv[3][3] = {
    {-0.5, 0.5, 0.5},
    { 0.5,-0.5, 0.5},
    { 0.5, 0.5,-0.5},
  };
  static double M_fcc_inv[3][3] = {
    { 0.0, 0.5, 0.5},
    { 0.5, 0.0, 0.5},
    { 0.5, 0.5, 0.0},
  };
  static double M_ac_inv[3][3] = {
    { 1.0, 0.0, 0.0},
    { 0.0, 0.5,-0.5},
    { 0.0, 0.5, 0.5},
  };
  static double M_bc_inv[3][3] = {
    { 0.5, 0.0,-0.5},
    { 0.0, 1.0, 0.0},
    { 0.5, 0.0, 0.5},
  };
  static double M_cc_inv[3][3] = {
    { 0.5, 0.5, 0.0},
    {-0.5, 0.5, 0.0},
    { 0.0, 0.0, 1.0},
  };
  static double M_rc_inv[3][3] = {
    { 2./3,-1./3,-1./3},
    { 1./3, 1./3,-2./3},
    { 1./3, 1./3, 1./3},
  };
  
  /* See: R. W. Grosse-Kunstleve, Acta Cryst. (1999). A55, 383-395 */
  /* */
  /* S is the Smith normal form, S = U*M*V */
  /* M = R - I, where R is virtical stack of generator rotation matrices */
  /* and I is virtical stack of identity matrices. */
  /* */
  /* Origin shift is given by (V*S+*U)*dw where dw is translations */
  /* corresponding to those obtained by symmetry finder and from */
  /* International table. */
  /* S+ is given by the operations of inversion of the diagonal elements */
  /* and then transpose. */
  
  static double tricli_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
  };
  
  static int tricli_generators[][3][9] =
  {
    { /* 1 */
      {  1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
  };
  
  static double monocli_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
  };
  
  static double monocli_A_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 1.0/2, 0,  0,  0,  0, },
    },
    { /* 9 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 1, 0, 0, -1.0/2, 0,  0,  0,  0, },
    },
  };
  
  static double monocli_B_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
  };
  
  static double monocli_C_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
  };
  
  static double monocli_I_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 1.0/2, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1, 1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      { -1, 0, 1.0/2,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, 1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 1.0/2, 0, -1.0/2,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1, -1.0/2, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, -1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
  };
  
  static int monocli_generators[][3][9] =
  {
    { /* 1 */
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  -1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  -1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  -1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 9 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
  };
  
  static double ortho_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { -1.0/2, 0, 0, 0, 0, 0, 0, 0, 0, },
      { 0, 0, 0, 0, -1.0/2, 0, 0, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static double ortho_F_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 1.0/2, 0, -1.0/4, -1.0/4, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, -1.0/4, -1.0/4, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, -1.0/4, 3.0/4, 0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 1, 0, 0, -1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, -1, 0, 0, 1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, -1, 0, 0, 1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 1, 0, 0, 0, 0, 0, -1.0/2, 0, },
      { 0, 0, 0, 0, 0, 0, 0, -1.0/2, 0, },
      { 0, 0, 0, 0, 1, 0, 0, -1.0/2, 0, },
    },
  };
  
  static double ortho_I_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1.0/2, -1.0/2, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, 0, 1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, 0, 0, 0, 1.0/2, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, -1.0/2, 0, -1, 0, 1.0/2, 0, 0, },
      { 0, 0, -1.0/2, 0, 0, 0, 0, 0, 0, },
    },
  };
  static double ortho_A_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 1, 0, 0, -1.0/2, 0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 1.0/2, 0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { -1.0/2, 0, 0, 0, 0, 0, 0, 0, 0, },
      { 0, -1, 0, 0, 0, -1.0/2, 0, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static double ortho_B_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, 0, -1.0/2, 0, 0, 0, 0, },
      { 0, 0, 0, -1, 0, 0, 1.0/2, 0, 0, },
    },
  };
  
  static double ortho_C_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static int ortho_generators[][3][9] =
  {
    { /* 1 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double tetra_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { -1.0/2, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1, 0, 0, 0, 1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 1.0/2, -1.0/2, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { -1, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 8 */
      { -1, 0, 0, 0, 1.0/2, 0, 0, 0, 0, },
      { 0, 0, 0, 0, -1.0/2, 0, 0, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static double tetra_I_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, -1, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { -1, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 0, -1, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -3.0/4, 1.0/4, 1.0/4,  0,  0,  0,  0,  0,  0, },
      { -1.0/4, -1.0/4, -1.0/4,  0,  0,  0,  0,  0,  0, },
      { -1.0/2, 1.0/2, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, -1, -1, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -3.0/4, 1.0/4, 0, 0, 1.0/4, 0,  0,  0,  0, },
      { -1.0/4, -1.0/4, 0, 0, -1.0/4, 0,  0,  0,  0, },
      { -1.0/2, 1.0/2, 0, 0, -1.0/2, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 2, 0, -1.0/2, 0, -1,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { -1, 0, 0, 0, -1, 0, 1.0/2, 0, -1.0/2, },
      { 0, 0, 0, 0, -1, 0, 1.0/2, 0, -1.0/2, },
      { 1, 0, 0, 0, 1, 0, -1, 0, 0, },
    },
  };
  
  static int tetra_generators[][3][9] =
  {
    { /* 1 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  0, 1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  0, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  0, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, 1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double trigo_VSpU[][3][9] =
  {
    { /* 1 */
      { -2.0/3, 1.0/3, 0,  0,  0,  0,  0,  0,  0, },
      { -1.0/3, -1.0/3, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 1.0/3, 0, -2.0/3, 0, 0,  0,  0,  0, },
      { 0, -1.0/3, 0, -1.0/3, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 3 */
      { 0, -1, 0, -2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 4 */
      { 0, -1, 0, -2, 0, 0,  0,  0,  0, },
      { 0, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, 1.0/3, 0, -2.0/3, 0, 0,  0,  0,  0, },
      { 0, -1.0/3, 0, -1.0/3, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 1, -1, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 1, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 1, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 10 */
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 11 */
      { -6, 3, 0, 4, 0, 0,  0,  0,  0, },
      { -3, 1, 0, 2, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 12 */
      { 0, 1, 0, -2, 0, 0, 1, 0, 0, },
      { 0, -1, 0, 1, 0, 0, -1, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
    { /* 13 */
      { 0, -1, 0, -2, 0, 0, 0, 0, 0, },
      { 0, -1, 0, -1, 0, 0, 0, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static int trigo_generators[][3][9] =
  {
    { /* 1 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 9 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 10 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 11 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 12 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
    { /* 13 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double rhombo_h_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 1,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/2, -1, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, -1.0/2, 1.0/2,  0,  0,  0,  0,  0,  0, },
      { 1.0/2, -1.0/2, -1.0/2,  0,  0,  0,  0,  0,  0, },
      { -1.0/2, 1.0/2, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, -1, 0, 0, 1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 1, 0, -1, -1.0/2, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1.0/2, -1, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 1, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, -1, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, -1, 0, 0, 1.0/2, 0, 0, },
    },
  };
  
  static int rhombo_h_generators[][3][9] =
  {
    { /* 1 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  0, 1, 0, -1, 1, 0, 0, 0, -1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  0, -1, 0, 1, -1, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double rhombo_p_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 1,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/2, -1, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 0, 0, -1, 0, 0,  0,  0,  0, },
      { 0, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/2, -1.0/2, 1.0/2,  0,  0,  0,  0,  0,  0, },
      { 1.0/2, -1.0/2, -1.0/2,  0,  0,  0,  0,  0,  0, },
      { -1.0/2, 1.0/2, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { 0, -1, 0, 0, 1.0/2, 0,  0,  0,  0, },
      { 0, 0, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, 1, 0, -1, -1.0/2, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -1.0/2, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1.0/2, -1, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/2, 0, 0, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 1, -1.0/2, 0, 0,  0,  0,  0, },
      { 1, 0, 0, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, -1, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, -1, 0, 0, 1.0/2, 0, 0, },
    },
  };
  
  static int rhombo_p_generators[][3][9] =
  {
    { /* 1 */
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  0, 0, -1, 0, -1, 0, -1, 0, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  0, 0, 1, 0, 1, 0, 1, 0, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {  0, 0, -1, 0, -1, 0, -1, 0, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {  0, 0, 1, 0, 1, 0, 1, 0, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  0, 0, -1, 0, -1, 0, -1, 0, 0, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double hexa_VSpU[][3][9] =
  {
    { /* 1 */
      { -1, 1, 0,  0,  0,  0,  0,  0,  0, },
      { -1, 0, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, 0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      { 1, 0, 0, -1, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 3 */
      { -1, 0, 0, -1, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { -1.0/3, -1.0/3, 0,  0,  0,  0,  0,  0,  0, },
      { 1.0/3, -2.0/3, 0,  0,  0,  0,  0,  0,  0, },
      { 0, 0, -1.0/2,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      { -1.0/3, 0, 0, -1.0/3, 0, 0,  0,  0,  0, },
      { 1.0/3, 0, 0, -2.0/3, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -1, 0, 0, 1, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { -1, 1, 0, 0, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 0, 0, 0, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 8 */
      { 1, 0, 0, -1, 0, 0, 0, 0, 0, },
      { -1, 0, 0, 0, 0, 0, 0, 0, 0, },
      { 0, 0, 0, 0, 0, -1.0/2, 0, 0, 0, },
    },
  };
  
  static int hexa_generators[][3][9] =
  {
    { /* 1 */
      {  1, -1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  1, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  1, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  -1, 1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  -1, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  -1, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  0, 1, 0, 1, 0, 0, 0, 0, 1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  1, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  1, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, -1, 0, -1, 0, 0, 0, 0, -1, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  static double cubic_VSpU[][3][9] =
  {
    { /* 1 */
      { -1.0/2, 1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, 1.0/2, 0, 1, 0, 0,  0,  0,  0, },
    },
    { /* 2 */
      { -1.0/2, 1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, -1.0/2, 0, -1, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, -1.0/2, 0, -1, 0, -1,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, -1,  0,  0,  0, },
    },
    { /* 4 */
      { 0, -1.0/2, 0, -1, 0, 1,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 0, 1.0/2, 0, 0, 0, -1,  0,  0,  0, },
    },
    { /* 5 */
      { -1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { -1.0/2, -1.0/2, 0, 1, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, -1.0/2, 0, 0, 0, 0,  0,  0,  0, },
      { 1.0/2, 1.0/2, 0, -1, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, 0, -1.0/2, -1, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 1,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 1.0/2, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 1.0/2, 0, 0, -1,  0,  0,  0, },
      { 0, 0, -1.0/2, 0, 0, 0,  0,  0,  0, },
    },
    { /* 9 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, -1, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 0, -1.0/2, 0, 0, },
    },
    { /* 10 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 1, -1.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 0, -1.0/2, 0, 0, },
    },
  };
  
  static double cubic_F_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, 0, -1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, 1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 2 */
      { 0, 1.0/2, -1.0/2, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1.0/2, 1.0/2, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1.0/2, -1.0/2, 0, 1.0/2, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 0, 1.0/4, -1.0/4, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -3.0/4, -1.0/4, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 1.0/4, -1.0/4, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 4 */
      { 0, 1.0/2, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, 0, -1.0/2, 0,  0,  0,  0, },
      { 0, -1.0/2, 0, -1, 1.0/2, 0,  0,  0,  0, },
    },
    { /* 5 */
      { -1.0/4, 1.0/4, 0, -1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/4, -3.0/4, 0, 1.0/2, 0, 0,  0,  0,  0, },
      { -1.0/4, 1.0/4, 0, 1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { 0, 0, 1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, -1, -1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
      { 0, 0, -1.0/2, -1.0/2, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { 0, -1.0/2, 0, -1.0/2, 0, -1.0/2,  0,  0,  0, },
      { 0, -1.0/2, 0, 1.0/2, 0, 1.0/2,  0,  0,  0, },
      { 0, -1.0/2, 0, 1.0/2, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 8 */
      { 0, 0, 1.0/2, -1.0/2, 0, 1.0/2,  0,  0,  0, },
      { 0, 0, 1.0/2, 1.0/2, 0, -1.0/2,  0,  0,  0, },
      { 0, 0, -1.0/2, -1.0/2, 0, -1.0/2,  0,  0,  0, },
    },
    { /* 9 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, -1, 0, -1, 0, 0, 1.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 0, -1.0/2, 0, 0, },
    },
    { /* 10 */
      { 0, 0, 0, 0, 0, 0, -1.0/2, 0, 0, },
      { 0, 0, -1, -2, 0, 0, 3.0/2, 0, 0, },
      { 0, 0, 0, 1, 0, 0, -1.0/2, 0, 0, },
    },
  };
  
  static double cubic_I_VSpU[][3][9] =
  {
    { /* 1 */
      { 0, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 1, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 0, -1, 0, 0, 0, 0,  0,  0,  0, },
    },
    { /* 2 */
      { -1, 0, 1, -1, 0, 0,  0,  0,  0, },
      { 0, 0, 1, -1, 0, 0,  0,  0,  0, },
      { 1, 0, -1, 0, 0, 0,  0,  0,  0, },
    },
    { /* 3 */
      { 1, 0, -1, -2, 0, -1,  0,  0,  0, },
      { 1, 0, -1, -1, 0, 0,  0,  0,  0, },
      { 1, 0, -1, -1, 0, -1,  0,  0,  0, },
    },
    { /* 4 */
      { -1, 0, 1, 0, 0, -1,  0,  0,  0, },
      { 3, 0, -1, -3, 0, 2,  0,  0,  0, },
      { -3, 0, 1, 3, 0, -3,  0,  0,  0, },
    },
    { /* 5 */
      { 1, 2, -5, -7, 0, 0,  0,  0,  0, },
      { 1, 1, -4, -5, 0, 0,  0,  0,  0, },
      { 0, 1, -2, -2, 0, 0,  0,  0,  0, },
    },
    { /* 6 */
      { -2, 1, 0, 1, 0, 0,  0,  0,  0, },
      { 1, -1, 0, -1, 0, 0,  0,  0,  0, },
      { 2, -1, 0, -2, 0, 0,  0,  0,  0, },
    },
    { /* 7 */
      { -1, 0, 0, 0, 0, -1,  0,  0,  0, },
      { -1, 0, 0, 1, 0, 0,  0,  0,  0, },
      { -1, 0, 0, 1, 0, -1,  0,  0,  0, },
    },
    { /* 8 */
      { 1, 0, 0, 0, -2, 1,  0,  0,  0, },
      { -1, 0, 0, 0, 1, -1,  0,  0,  0, },
      { 1, 0, 0, 0, -1, 0,  0,  0,  0, },
    },
    { /* 9 */
      { -2, -3, 0, -1, 0, 0, 1, -1, 1, },
      { -1, -3, 0, -1, 0, 0, 1, -1, 1, },
      { 0, -1, 0, 0, 0, 0, 0, 0, 0, },
    },
    { /* 10 */
      { -1, 0, 0, 0, 0, -1, 0, -1, 1, },
      { -1, 0, 0, 1, 0, 0, 0, -1, 1, },
      { -1, 0, 0, 1, 0, -1, 0, -1, 1, },
    },
  };
  
  static int cubic_generators[][3][9] =
  {
    { /* 1 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 2 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 3 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 4 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 5 */
      {  0, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 6 */
      {  0, 1, 0, -1, 0, 0, 0, 0, -1, },
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 7 */
      {  1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 8 */
      {  1, 0, 0, 0, 1, 0, 0, 0, -1, },
      {  0, 0, -1, -1, 0, 0, 0, -1, 0, },
      {   0,  0,  0,  0,  0,  0,  0,  0,  0, },
    },
    { /* 9 */
      {  0, -1, 0, 1, 0, 0, 0, 0, 1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
    { /* 10 */
      {  -1, 0, 0, 0, -1, 0, 0, 0, 1, },
      {  0, 0, 1, 1, 0, 0, 0, 1, 0, },
      {  -1, 0, 0, 0, -1, 0, 0, 0, -1, },
    },
  };
  
  
  static int find_hall_symbol(double origin_shift[3],
                              SPGCONST double bravais_lattice[3][3],
                              const int hall_number,
                              const Centering centering,
                              const Symmetry *symmetry,
                              const double symprec);
  static int is_hall_symbol(double shift[3],
                            const int hall_number,
                            SPGCONST double primitive_lattice[3][3],
                            const Symmetry *symmetry,
                            const Centering centering,
                            SPGCONST int generators[3][9],
                            SPGCONST double VSpU[3][9],
                            const double symprec);
  static int is_hall_symbol_cubic(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec);
  static int is_hall_symbol_hexa(double shift[3],
                                 const int hall_number,
                                 SPGCONST double primitive_lattice[3][3],
                                 const Symmetry *symmetry,
                                 const double symprec);
  static int is_hall_symbol_rhombo(double shift[3],
                                   const int hall_number,
                                   SPGCONST double primitive_lattice[3][3],
                                   const Symmetry *symmetry,
                                   const double symprec);
  static int is_hall_symbol_trigonal(double shift[3],
                                     const int hall_number,
                                     SPGCONST double primitive_lattice[3][3],
                                     const Symmetry *symmetry,
                                     const double symprec);
  static int is_hall_symbol_tetra(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec);
  static int is_hall_symbol_ortho(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec);
  static int is_hall_symbol_monocli(double shift[3],
                                    const int hall_number,
                                    SPGCONST double primitive_lattice[3][3],
                                    const Symmetry *symmetry,
                                    const Centering centering,
                                    const double symprec);
  static int is_hall_symbol_tricli(double shift[3],
                                   const int hall_number,
                                   SPGCONST double primitive_lattice[3][3],
                                   const Symmetry *symmetry,
                                   const double symprec);
  static int get_translations(double trans[3][3],
                              const Symmetry *symmetry,
                              SPGCONST int rot[3][3][3]);
  static void transform_translation(double trans_reduced[3],
                                    const Centering centering,
                                    const double trans[3]);
  static void transform_rotation(double rot_reduced[3][3],
                                 const Centering centering,
                                 SPGCONST int rot[3][3]);
  static int get_origin_shift(double shift[3],
                              const int hall_number,
                              SPGCONST int rot[3][3][3],
                              SPGCONST double trans[3][3],
                              const Centering centering,
                              SPGCONST double VSpU[3][9]);
  static void unpack_generators(int rot[3][3][3], int generators[3][9]);
  static int set_dw(double dw[3],
                    const int operation_index[2],
                    SPGCONST int rot[3][3],
                    const double trans[3],
                    const Centering centering);
  static int is_match_database(const int hall_number,
                               const double shift[3],
                               SPGCONST double primitive_lattice[3][3],
                               const Centering centering,
                               const Symmetry *symmetry,
                               const double symprec);
  
  
  int hal_match_hall_symbol_db(double origin_shift[3],
                               SPGCONST double bravais_lattice[3][3],
                               const int hall_number,
                               const Centering centering,
                               const Symmetry *symmetry,
                               const double symprec)
  {
    return find_hall_symbol(origin_shift,
                            bravais_lattice,
                            hall_number,
                            centering,
                            symmetry,
                            symprec);
  }
  
  static int find_hall_symbol(double origin_shift[3],
                              SPGCONST double bravais_lattice[3][3],
                              const int hall_number,
                              const Centering centering,
                              const Symmetry *symmetry,
                              const double symprec)
  {
    double primitive_lattice[3][3];
  
    switch (centering) {
    case PRIMITIVE:
      mat_copy_matrix_d3(primitive_lattice, bravais_lattice);
      break;
    case BODY:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_bcc_inv);
      break;
    case FACE:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_fcc_inv);
      break;
    case A_FACE:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_ac_inv);
      break;
    case B_FACE:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_bc_inv);
      break;
    case C_FACE:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_cc_inv);
      break;
    case R_CENTER:
      mat_multiply_matrix_d3(primitive_lattice, bravais_lattice, M_rc_inv);
      break;
    default:
      break;
    }
  
    /* CUBIC IT: 195-230, Hall: 489-530 */
    if (489 <= hall_number && hall_number <= 530) {
      if (is_hall_symbol_cubic(origin_shift,
                               hall_number,
                               primitive_lattice,
                               symmetry,
                               centering,
                               symprec)) {goto found;}
      return 0;
    }
  
    /* HEXA, IT: 168-194, Hall: 462-488 */
    if (462 <= hall_number && hall_number <= 488) {
      if (is_hall_symbol_hexa(origin_shift,
                              hall_number,
                              primitive_lattice,
                              symmetry,
                              symprec)) {goto found;}
      return 0;
    }
  
    /* TRIGO, IT: 143-167, Hall: 430-461 */
    if (430 <= hall_number && hall_number <= 461) {
      if (hall_number == 433 ||
          hall_number == 434 ||
          hall_number == 436 ||
          hall_number == 437 ||
          hall_number == 444 ||
          hall_number == 445 ||
          hall_number == 450 ||
          hall_number == 451 ||
          hall_number == 452 ||
          hall_number == 453 ||
          hall_number == 458 ||
          hall_number == 459 ||
          hall_number == 460 ||
          hall_number == 461) {
        if (is_hall_symbol_rhombo(origin_shift,
                                  hall_number,
                                  primitive_lattice,
                                  symmetry,
                                  symprec)) {goto found;}
      } else {
        if (is_hall_symbol_trigonal(origin_shift,
                                    hall_number,
                                    primitive_lattice,
                                    symmetry,
                                    symprec)) {goto found;}
      }
      return 0;
    }
  
    /* TETRA, IT: 75-142, Hall: 349-429 */
    if (349 <= hall_number && hall_number <= 429) {
      if (is_hall_symbol_tetra(origin_shift,
                               hall_number,
                               primitive_lattice,
                               symmetry,
                               centering,
                               symprec)) {goto found;}
      return 0;
    }
  
    /* ORTHO, IT: 16-74, Hall: 108-348 */
    if (108 <= hall_number && hall_number <= 348) {
      if (is_hall_symbol_ortho(origin_shift,
                               hall_number,
                               primitive_lattice,
                               symmetry,
                               centering,
                               symprec)) {goto found;}
      return 0;
    }
  
    /* MONOCLI, IT: 3-15, Hall: 3-107 */
    if (3 <= hall_number && hall_number <= 107) {
      if (is_hall_symbol_monocli(origin_shift,
                                 hall_number,
                                 primitive_lattice,
                                 symmetry,
                                 centering,
                                 symprec)) {goto found;}
      return 0;
    }
  
    /* TRICLI, IT: 1-2, Hall: 1-2 */
    if (1 <= hall_number && hall_number <= 2) {
      if (is_hall_symbol_tricli(origin_shift,
                                hall_number,
                                primitive_lattice,
                                symmetry,
                                symprec)) {goto found;}
      return 0;
    }
  
    return 0;
  
   found:
    switch (centering) {
    case PRIMITIVE:
      break;
    case BODY:
      mat_multiply_matrix_vector_d3(origin_shift, M_bcc_inv, origin_shift);
      break;
    case FACE:
      mat_multiply_matrix_vector_d3(origin_shift, M_fcc_inv, origin_shift);
      break;
    case A_FACE:
      mat_multiply_matrix_vector_d3(origin_shift, M_ac_inv, origin_shift);
      break;
    case B_FACE:
      mat_multiply_matrix_vector_d3(origin_shift, M_bc_inv, origin_shift);
      break;
    case C_FACE:
      mat_multiply_matrix_vector_d3(origin_shift, M_cc_inv, origin_shift);
      break;
    case R_CENTER:
      mat_multiply_matrix_vector_d3(origin_shift, M_rc_inv, origin_shift);
      break;
    default:
      break;
    }
    return 1;
  }
  
  static int is_hall_symbol_cubic(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec)
  {
    int i;
  
    for (i = 0; i < 10; i++) {
      if (centering == PRIMITIVE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           cubic_generators[i],
                           cubic_VSpU[i],
                           symprec)) {goto found;}
      }
  
      if (centering == BODY) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           cubic_generators[i],
                           cubic_I_VSpU[i],
                           symprec)) {goto found;}
      }
  
      if (centering == FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           cubic_generators[i],
                           cubic_F_VSpU[i],
                           symprec)) {goto found;}
      }
    }
  
    return 0;
  
  found:
    return 1;
  }
  
  static int is_hall_symbol_hexa(double shift[3],
                                 const int hall_number,
                                 SPGCONST double primitive_lattice[3][3],
                                 const Symmetry *symmetry,
                                 const double symprec)
  {
    int i;
  
    for (i = 0; i < 8; i++) {
      if (is_hall_symbol(shift,
                         hall_number,
                         primitive_lattice,
                         symmetry,
                         PRIMITIVE,
                         hexa_generators[i],
                         hexa_VSpU[i],
                         symprec)) {return 1;}
    }
  
    return 0;
  }
  
  static int is_hall_symbol_trigonal(double shift[3],
                                     const int hall_number,
                                     SPGCONST double primitive_lattice[3][3],
                                     const Symmetry *symmetry,
                                     const double symprec)
  {
    int i;
  
    for (i = 0; i < 13; i++) {
      if (is_hall_symbol(shift,
                         hall_number,
                         primitive_lattice,
                         symmetry,
                         PRIMITIVE,
                         trigo_generators[i],
                         trigo_VSpU[i],
                         symprec)) {return 1;}
    }
  
    return 0;
  }
  
  static int is_hall_symbol_rhombo(double shift[3],
                                   const int hall_number,
                                   SPGCONST double primitive_lattice[3][3],
                                   const Symmetry *symmetry,
                                   const double symprec)
  {
    int i;
  
    if (hall_number == 433 ||
        hall_number == 436 ||
        hall_number == 444 ||
        hall_number == 450 ||
        hall_number == 452 ||
        hall_number == 458 ||
        hall_number == 460) {
      /* hP */
      for (i = 0; i < 8; i++) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           R_CENTER,
                           rhombo_h_generators[i],
                           rhombo_h_VSpU[i],
                           symprec)) {return 1;}
      }
    } else {
      /* hR */
      for (i = 0; i < 8; i++) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           PRIMITIVE,
                           rhombo_p_generators[i],
                           rhombo_p_VSpU[i],
                           symprec)) {return 1;}
      }
    }
  
    return 0;
  }
  
  static int is_hall_symbol_tetra(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec)
  {
    int i;
  
    for (i = 0; i < 8; i++) {
      if (centering==PRIMITIVE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           tetra_generators[i],
                           tetra_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering==BODY) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           tetra_generators[i],
                           tetra_I_VSpU[i],
                           symprec)) {return 1;}
      }
    }
  
    return 0;
  }
  
  static int is_hall_symbol_ortho(double shift[3],
                                  const int hall_number,
                                  SPGCONST double primitive_lattice[3][3],
                                  const Symmetry *symmetry,
                                  const Centering centering,
                                  const double symprec)
  {
    int i;
  
    for (i = 0; i < 5; i++) {
      if (centering == PRIMITIVE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == BODY) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_I_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_F_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == A_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_A_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == B_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_B_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == C_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           ortho_generators[i],
                           ortho_C_VSpU[i],
                           symprec)) {return 1;}
      }
    }
  
    return 0;
  }
  
  static int is_hall_symbol_monocli(double shift[3],
                                    const int hall_number,
                                    SPGCONST double primitive_lattice[3][3],
                                    const Symmetry *symmetry,
                                    const Centering centering,
                                    const double symprec)
  {
    int i;
  
    for (i = 0; i < 9; i++) {
      if (centering == PRIMITIVE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           monocli_generators[i],
                           monocli_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == A_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           monocli_generators[i],
                           monocli_A_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == B_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           monocli_generators[i],
                           monocli_B_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == C_FACE) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           monocli_generators[i],
                           monocli_C_VSpU[i],
                           symprec)) {return 1;}
      }
  
      if (centering == BODY) {
        if (is_hall_symbol(shift,
                           hall_number,
                           primitive_lattice,
                           symmetry,
                           centering,
                           monocli_generators[i],
                           monocli_I_VSpU[i],
                           symprec)) {return 1;}
      }
    }
  
    return 0;
  }
  
  static int is_hall_symbol_tricli(double shift[3],
                                   const int hall_number,
                                   SPGCONST double primitive_lattice[3][3],
                                   const Symmetry *symmetry,
                                   const double symprec)
  {
    int i;
  
    for (i = 0; i < 2; i++) {
      if (is_hall_symbol(shift,
                         hall_number,
                         primitive_lattice,
                         symmetry,
                         PRIMITIVE,
                         tricli_generators[i],
                         tricli_VSpU[i],
                         symprec)) {return 1;}
    }
  
    return 0;
  }
  
  static void unpack_generators(int rot[3][3][3], int generators[3][9])
  {
    int i, j, k;
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++) {
        for (k = 0; k < 3; k++) {
          rot[i][j][k] = generators[i][j*3+k];
        }
      }
    }
  }
  
  static int is_hall_symbol(double shift[3],
                            const int hall_number,
                            SPGCONST double primitive_lattice[3][3],
                            const Symmetry *symmetry,
                            const Centering centering,
                            SPGCONST int generators[3][9],
                            SPGCONST double VSpU[3][9],
                            const double symprec)
  {
    int is_origin_shift;
    int operation_index[2];
    int rot[3][3][3];
    double trans[3][3];
  
    debug_print("[line %d, %s]\n", __LINE__, __FILE__);
    debug_print("primitive lattice\n");
    debug_print_matrix_d3(primitive_lattice);
  
    spgdb_get_operation_index(operation_index, hall_number);
  
    if (! (operation_index[0] == symmetry->size)) {
      goto not_found;
    }
  
    unpack_generators(rot, generators);
    if (get_translations(trans, symmetry, rot)) {
      is_origin_shift = get_origin_shift(shift,
                                         hall_number,
                                         rot,
                                         trans,
                                         centering,
                                         VSpU);
  
      if (is_origin_shift) {
        if (is_match_database(hall_number,
                              shift,
                              primitive_lattice,
                              centering,
                              symmetry,
                              symprec)) {goto found;}
      }
    } else {
      goto not_found;
    }
  
   not_found:
    return 0;
  
   found:
    debug_print("[line %d, %s]\n", __LINE__, __FILE__);
    debug_print("origin shift\n");
    debug_print_vector_d3(shift);
  
    return 1;
  }
  
  static int get_translations(double trans[3][3],
                              const Symmetry *symmetry,
                              SPGCONST int rot[3][3][3])
  {
    int i, j;
    int is_found;
    static SPGCONST int zero[3][3] = { { 0, 0, 0 },
                                       { 0, 0, 0 },
                                       { 0, 0, 0 }, };
  
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++) {
        trans[i][j] = 0;
      }
    }
  
    for (i = 0; i < 3; i++) {
      if (mat_check_identity_matrix_i3(rot[i], zero)) { continue; }
      is_found = 0;
      for (j = 0; j < symmetry->size; j++) {
        if (mat_check_identity_matrix_i3(symmetry->rot[j], rot[i])) {
          mat_copy_vector_d3(trans[i], symmetry->trans[j]);
          is_found = 1;
          break;
        }
      }
      if (! is_found) {
        goto not_found;
      }
    }
  
    return 1;
  
   not_found:
  
    return 0;
  }
  
  static void transform_translation(double trans_reduced[3],
                                    const Centering centering,
                                    const double trans[3])
  {
    int i;
  
    switch (centering) {
    case PRIMITIVE:
      mat_copy_vector_d3(trans_reduced, trans);
      break;
    case BODY:
      mat_multiply_matrix_vector_id3(trans_reduced, M_bcc, trans);
      break;
    case FACE:
      mat_multiply_matrix_vector_id3(trans_reduced, M_fcc, trans);
      break;
    case A_FACE:
      mat_multiply_matrix_vector_id3(trans_reduced, M_ac, trans);
      break;
    case B_FACE:
      mat_multiply_matrix_vector_id3(trans_reduced, M_bc, trans);
      break;
    case C_FACE:
      mat_multiply_matrix_vector_id3(trans_reduced, M_cc, trans);
      break;
    case R_CENTER:
      mat_multiply_matrix_vector_id3(trans_reduced, M_rc, trans);
      break;
    default:
      break;
    }
  
    for (i = 0; i < 3; i++) {
      trans_reduced[i] = mat_Dmod1(trans_reduced[i]);
    }
  }
  
  static void transform_rotation(double rot_reduced[3][3],
                                 const Centering centering,
                                 SPGCONST int rot[3][3])
  {
    mat_cast_matrix_3i_to_3d(rot_reduced, rot);
    if (centering != PRIMITIVE) {
      switch (centering) {
      case BODY:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_bcc_inv, 0);
        break;
      case FACE:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_fcc_inv, 0);
        break;
      case A_FACE:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_ac_inv, 0);
        break;
      case B_FACE:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_bc_inv, 0);
        break;
      case C_FACE:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_cc_inv, 0);
        break;
      case R_CENTER:
        mat_get_similar_matrix_d3(rot_reduced, rot_reduced, M_rc_inv, 0);
        break;
      default:
        break;
      }
    }
  }
  
  static int get_origin_shift(double shift[3],
                              const int hall_number,
                              SPGCONST int rot[3][3][3],
                              SPGCONST double trans[3][3],
                              const Centering centering,
                              SPGCONST double VSpU[3][9])
  {
    int i, j;
    int operation_index[2];
    double dw[9], tmp_dw[3];
  
    spgdb_get_operation_index(operation_index, hall_number);
  
    /* The obtained dw is reduced to that of primitve cell by centerings. */
    for (i = 0; i < 3; i++) {
      /* Zero matrix is the sign to set dw 0 */
      if (mat_get_determinant_i3(rot[i]) == 0) {
        for (j = 0; j < 3; j++) {
          dw[i * 3 + j] = 0;
        }
      } else {
        if (set_dw(tmp_dw, operation_index, rot[i], trans[i], centering)) {
          for (j = 0; j < 3; j++) {
            dw[i * 3 + j] = tmp_dw[j];
          }
        } else {
          goto not_found;
        }
      }
    }
  
    /* VSpU*dw is given for the primitive cell if there is centering. */
    for (i = 0; i < 3; i++) {
      shift[i] = 0;
      for (j = 0; j < 9; j++) {
        shift[i] += VSpU[i][j] * dw[j];
      }
      shift[i] = mat_Dmod1(shift[i]);
    }
  
    return 1;
  
   not_found:
    return 0;
  }
  
  static int set_dw(double dw[3],
                    const int operation_index[2],
                    SPGCONST int rot[3][3],
                    const double trans[3],
                    const Centering centering)
  {
    int i, j;
    int rot_db[3][3];
    double trans_db[3], trans_prim[3], trans_db_prim[3];
  
    transform_translation(trans_prim, centering, trans);
    for (i = 0; i < operation_index[0]; i++) {
      spgdb_get_operation(rot_db, trans_db, operation_index[1] + i);
      transform_translation(trans_db_prim, centering, trans_db);
      if (mat_check_identity_matrix_i3(rot_db, rot)) {
        for (j = 0; j < 3; j++) {
          dw[j] = trans_prim[j] - trans_db_prim[j];
          dw[j] = mat_Dmod1(dw[j]);
        }
        goto found;
      }
    }
  
    /* Not found */
    return 0;
  
   found:
    return 1;
  }
  
  static int is_match_database(const int hall_number,
                               const double origin_shift[3],
                               SPGCONST double primitive_lattice[3][3],
                               const Centering centering,
                               const Symmetry *symmetry,
                               const double symprec)
  {
    int i, j, k, is_found;
    int operation_index[2];
    int rot_db[3][3];
    int found_list[192];
    double trans_db[3], trans_db_prim[3], trans_prim[3], diff[3], shift_rot[3];
    double rot_prim[3][3];
  
    spgdb_get_operation_index(operation_index, hall_number);
  
    for (i = 0; i < symmetry->size; i++) {found_list[i] = 0;}
    for (i = 0; i < symmetry->size; i++) {
      is_found = 0;
      for (j = 0; j < operation_index[0]; j++) {
        spgdb_get_operation(rot_db, trans_db, operation_index[1] + j);
        if (mat_check_identity_matrix_i3(symmetry->rot[i], rot_db)) {
          transform_translation(trans_db_prim, centering, trans_db);
          transform_translation(trans_prim, centering, symmetry->trans[i]);
          transform_rotation(rot_prim, centering, rot_db);
          for (k = 0; k < 3; k++) {
            diff[k] = trans_prim[k] - trans_db_prim[k] + origin_shift[k];
          }
          mat_multiply_matrix_vector_d3(shift_rot, rot_prim, origin_shift);
          if (cel_is_overlap(diff, shift_rot, primitive_lattice, symprec)) {
            if (! found_list[j]) {
              found_list[j] = 1;
              is_found = 1;
              break;
            }
          }
        }
      }
  
      if (! is_found) {
        goto not_found;
      }
    }
  
    return 1;
  
   not_found:
    return 0;
  }
  