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