hall_symbol.c

/* 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 phonopy 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 "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,
			    SPGCONST Symmetry *symmetry,
			    const double symprec);
static int is_hall_symbol(double shift[3],
			  const int hall_number,
			  SPGCONST double primitive_lattice[3][3],
			  SPGCONST Symmetry *symmetry,
			  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],
				SPGCONST 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],
			       SPGCONST Symmetry *symmetry,
			       const double symprec);
static int is_hall_symbol_rhombo(double shift[3],
				 const int hall_number,
				 SPGCONST double primitive_lattice[3][3],
				 SPGCONST Symmetry *symmetry,
				 const double symprec);
static int is_hall_symbol_trigonal(double shift[3],
				   const int hall_number,
				   SPGCONST double primitive_lattice[3][3],
				   SPGCONST Symmetry *symmetry,
				   const double symprec);
static int is_hall_symbol_tetra(double shift[3],
				const int hall_number,
				SPGCONST double primitive_lattice[3][3],
				SPGCONST 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],
				SPGCONST 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],
				  SPGCONST 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],
				 SPGCONST Symmetry *symmetry,
				 const double symprec);
static int get_translations(double trans[3][3],
			    SPGCONST 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,
			     SPGCONST 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,
			     SPGCONST 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,
			    SPGCONST 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 (centering == R_CENTER) {
      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],
				SPGCONST 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)) {return 1;}
    }

    if (centering == BODY) {
      if (is_hall_symbol(shift,
			 hall_number,
			 primitive_lattice,
			 symmetry,
			 centering,
			 cubic_generators[i],
			 cubic_I_VSpU[i],
			 symprec)) {return 1;}
    }

    if (centering == FACE) {
      if (is_hall_symbol(shift,
			 hall_number,
			 primitive_lattice,
			 symmetry,
			 centering,
			 cubic_generators[i],
			 cubic_F_VSpU[i],
			 symprec)) {return 1;}
    }
  }

  return 0;
}

static int is_hall_symbol_hexa(double shift[3],
			       const int hall_number,
			       SPGCONST double primitive_lattice[3][3],
			       SPGCONST 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],
				   SPGCONST 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],
				 SPGCONST 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],
				SPGCONST 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],
				SPGCONST 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],
				  SPGCONST 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],
				 SPGCONST 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],
			  SPGCONST Symmetry *symmetry,
			  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];

  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:
  return 1;
}

static int get_translations(double trans[3][3],
			    SPGCONST 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,
			     SPGCONST 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;
}