text: /* Copyright (C) 2008 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 <string.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include "pointgroup.h"
  #include "symmetry.h"
  #include "mathfunc.h"
  
  #include "debug.h"
  
  #define NUM_ROT_AXES 73
  
  typedef struct {
    int table[10];
    char symbol[6];
    char schoenflies[4];
    Holohedry holohedry;
    Laue laue;
  } PointgroupType;
  
  static PointgroupType pointgroup_data[33] = {
    { /* 0 */
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      "     ",
      "   ",
      HOLOHEDRY_NONE,
      LAUE_NONE,
    },
    { /* 1 */
      {0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
      "1    ",
      "C1 ",
      TRICLI,
      LAUE1,
    },
    { /* 2 */
      {0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
      "-1   ",
      "Ci ",
      TRICLI,
      LAUE1,
    },
    { /* 3 */
      {0, 0, 0, 0, 0, 1, 1, 0, 0, 0},
      "2    ",
      "C2 ",
      MONOCLI,
      LAUE2M,
    },
    { /* 4 */
      {0, 0, 0, 1, 0, 1, 0, 0, 0, 0},
      "m    ",
      "Cs ",
      MONOCLI,
      LAUE2M,
    },
    { /* 5 */
      {0, 0, 0, 1, 1, 1, 1, 0, 0, 0},
      "2/m  ",
      "C2h",
      MONOCLI,
      LAUE2M,
    },
    { /* 6 */
      {0, 0, 0, 0, 0, 1, 3, 0, 0, 0},
      "222  ",
      "D2 ",
      ORTHO,
      LAUEMMM,
    },
    { /* 7 */
      {0, 0, 0, 2, 0, 1, 1, 0, 0, 0},
      "mm2  ",
      "C2v",
      ORTHO,
      LAUEMMM,
    },
    { /* 8 */
      {0, 0, 0, 3, 1, 1, 3, 0, 0, 0},
      "mmm  ",
      "D2h",
      ORTHO,
      LAUEMMM,
    },
    { /* 9 */
      {0, 0, 0, 0, 0, 1, 1, 0, 2, 0},
      "4    ",
      "C4 ",
      TETRA,
      LAUE4M,
    },
    { /* 10 */
      {0, 2, 0, 0, 0, 1, 1, 0, 0, 0},
      "-4   ",
      "S4 ",
      TETRA,
      LAUE4M,
    },
    { /* 11 */
      {0, 2, 0, 1, 1, 1, 1, 0, 2, 0},
      "4/m  ",
      "C4h",
      TETRA,
      LAUE4M,
    },
    { /* 12 */
      {0, 0, 0, 0, 0, 1, 5, 0, 2, 0},
      "422  ",
      "D4 ",
      TETRA,
      LAUE4MMM,
    },
    { /* 13 */
      {0, 0, 0, 4, 0, 1, 1, 0, 2, 0},
      "4mm  ",
      "C4v",
      TETRA,
      LAUE4MMM,
    },
    { /* 14 */
      {0, 2, 0, 2, 0, 1, 3, 0, 0, 0},
      "-42m ",
      "D2d",
      TETRA,
      LAUE4MMM,
    },
    { /* 15 */
      {0, 2, 0, 5, 1, 1, 5, 0, 2, 0},
      "4/mmm",
      "D4h",
      TETRA,
      LAUE4MMM,
    },
    { /* 16 */
      {0, 0, 0, 0, 0, 1, 0, 2, 0, 0},
      "3    ",
      "C3 ",
      TRIGO,
      LAUE3,
    },
    { /* 17 */
      {0, 0, 2, 0, 1, 1, 0, 2, 0, 0},
      "-3   ",
      "C3i",
      TRIGO,
      LAUE3,
    },
    { /* 18 */
      {0, 0, 0, 0, 0, 1, 3, 2, 0, 0},
      "32   ",
      "D3 ",
      TRIGO,
      LAUE3M,
    },
    { /* 19 */
      {0, 0, 0, 3, 0, 1, 0, 2, 0, 0},
      "3m   ",
      "C3v",
      TRIGO,
      LAUE3M,
    },
    { /* 20 */
      {0, 0, 2, 3, 1, 1, 3, 2, 0, 0},
      "-3m  ",
      "D3d",
      TRIGO,
      LAUE3M,
    },
    { /* 21 */
      {0, 0, 0, 0, 0, 1, 1, 2, 0, 2},
      "6    ",
      "C6 ",
      HEXA,
      LAUE6M,
    },
    { /* 22 */
      {2, 0, 0, 1, 0, 1, 0, 2, 0, 0},
      "-6   ",
      "C3h",
      HEXA,
      LAUE6M,
    },
    { /* 23 */
      {2, 0, 2, 1, 1, 1, 1, 2, 0, 2},
      "6/m  ",
      "C6h",
      HEXA,
      LAUE6M,
    },
    { /* 24 */
      {0, 0, 0, 0, 0, 1, 7, 2, 0, 2},
      "622  ",
      "D6 ",
      HEXA,
      LAUE6MMM,
    },
    { /* 25 */
      {0, 0, 0, 6, 0, 1, 1, 2, 0, 2},
      "6mm  ",
      "C6v",
      HEXA,
      LAUE6MMM,
    },
    { /* 26 */
      {2, 0, 0, 4, 0, 1, 3, 2, 0, 0},
      "-6m2 ",
      "D3h",
      HEXA,
      LAUE6MMM,
    },
    { /* 27 */
      {2, 0, 2, 7, 1, 1, 7, 2, 0, 2},
      "6/mmm",
      "D6h",
      HEXA,
      LAUE6MMM,
    },
    { /* 28 */
      {0, 0, 0, 0, 0, 1, 3, 8, 0, 0},
      "23   ",
      "T  ",
      CUBIC,
      LAUEM3,
    },
    { /* 29 */
      {0, 0, 8, 3, 1, 1, 3, 8, 0, 0},
      "m-3  ",
      "Th ",
      CUBIC,
      LAUEM3,
    },
    { /* 30 */
      {0, 0, 0, 0, 0, 1, 9, 8, 6, 0},
      "432  ",
      "O  ",
      CUBIC,
      LAUEM3M,
    },
    { /* 31 */
      {0, 6, 0, 6, 0, 1, 3, 8, 0, 0},
      "-43m ",
      "Td ",
      CUBIC,
      LAUEM3M,
    },
    { /* 32 */
      {0, 6, 8, 9, 1, 1, 9, 8, 6, 0},
      "m-3m ",
      "Oh ",
      CUBIC,
      LAUEM3M,
    }
  };
  
  static int identity[3][3] = {
    { 1, 0, 0},
    { 0, 1, 0},
    { 0, 0, 1},
  };
  
  static int inversion[3][3] = {
    {-1, 0, 0},
    { 0,-1, 0},
    { 0, 0,-1},
  };
  
  static int rot_axes[][3] = {
    { 1, 0, 0},
    { 0, 1, 0},
    { 0, 0, 1},
    { 0, 1, 1},
    { 1, 0, 1},
    { 1, 1, 0},
    { 0, 1,-1},
    {-1, 0, 1},
    { 1,-1, 0},
    { 1, 1, 1}, /* 10 */
    {-1, 1, 1},
    { 1,-1, 1},
    { 1, 1,-1},
    { 0, 1, 2},
    { 2, 0, 1},
    { 1, 2, 0},
    { 0, 2, 1},
    { 1, 0, 2},
    { 2, 1, 0},
    { 0,-1, 2}, /* 20 */
    { 2, 0,-1},
    {-1, 2, 0},
    { 0,-2, 1},
    { 1, 0,-2},
    {-2, 1, 0},
    { 2, 1, 1},
    { 1, 2, 1},
    { 1, 1, 2},
    { 2,-1,-1},
    {-1, 2,-1}, /* 30 */
    {-1,-1, 2},
    { 2, 1,-1},
    {-1, 2, 1},
    { 1,-1, 2},
    { 2,-1, 1}, /* 35 */
    { 1, 2,-1},
    {-1, 1, 2},
    { 3, 1, 2},
    { 2, 3, 1},
    { 1, 2, 3}, /* 40 */
    { 3, 2, 1},
    { 1, 3, 2},
    { 2, 1, 3},
    { 3,-1, 2},
    { 2, 3,-1}, /* 45 */
    {-1, 2, 3},
    { 3,-2, 1},
    { 1, 3,-2},
    {-2, 1, 3},
    { 3,-1,-2}, /* 50 */
    {-2, 3,-1},
    {-1,-2, 3},
    { 3,-2,-1},
    {-1, 3,-2},
    {-2,-1, 3}, /* 55 */
    { 3, 1,-2},
    {-2, 3, 1},
    { 1,-2, 3},
    { 3, 2,-1},
    {-1, 3, 2}, /* 60 */
    { 2,-1, 3},
    { 1, 1, 3},
    {-1, 1, 3},
    { 1,-1, 3},
    {-1,-1, 3}, /* 65 */
    { 1, 3, 1},
    {-1, 3, 1},
    { 1, 3,-1},
    {-1, 3,-1},
    { 3, 1, 1}, /* 70 */
    { 3, 1,-1},
    { 3,-1, 1},
    { 3,-1,-1},
  };
  
  static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
                                                const int num_rotations);
  static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym);
  static int get_pointgroup_class_table(int table[10],
                                        SPGCONST PointSymmetry * pointsym);
  static int get_rotation_type(SPGCONST int rot[3][3]);
  static int get_rotation_axis(SPGCONST int rot[3][3]);
  static int get_orthogonal_axis(int ortho_axes[],
                                 SPGCONST int proper_rot[3][3],
                                 const int rot_order);
  static int laue2m(int axes[3],
                    SPGCONST PointSymmetry * pointsym);
  
  #ifdef SPGDEBUG
  static int lauemmm(int axes[3],
                     SPGCONST PointSymmetry * pointsym);
  static int laue4m(int axes[3],
                    SPGCONST PointSymmetry * pointsym);
  static int laue4mmm(int axes[3],
                      SPGCONST PointSymmetry * pointsym);
  static int laue3(int axes[3],
                   SPGCONST PointSymmetry * pointsym);
  static int laue3m(int axes[3],
                    SPGCONST PointSymmetry * pointsym);
  static int lauem3m(int axes[3],
                     SPGCONST PointSymmetry * pointsym);
  #endif
  
  static int laue_one_axis(int axes[3],
                           SPGCONST PointSymmetry * pointsym,
                           const int rot_order);
  static int lauennn(int axes[3],
                     SPGCONST PointSymmetry * pointsym,
                     const int rot_order);
  static int get_axes(int axes[3],
                      const Laue laue,
                      SPGCONST PointSymmetry * pointsym);
  static void get_proper_rotation(int prop_rot[3][3],
                                  SPGCONST int rot[3][3]);
  static void set_transformation_matrix(int tmat[3][3],
                                        const int axes[3]);
  static int is_exist_axis(const int axis_vec[3], const int axis_index);
  static void sort_axes(int axes[3]);
  
  /* Retrun pointgroup.number = 0 if failed */
  Pointgroup ptg_get_transformation_matrix(int transform_mat[3][3],
                                           SPGCONST int rotations[][3][3],
                                           const int num_rotations)
  {
    int i, j, pg_num;
    int axes[3];
    PointSymmetry pointsym;
    Pointgroup pointgroup;
  
    debug_print("ptg_get_transformation_matrix:\n");
  
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++) {
        transform_mat[i][j] = 0;
      }
    }
  
    pg_num = get_pointgroup_number_by_rotations(rotations, num_rotations);
  
    if (pg_num > 0) {
      pointgroup = ptg_get_pointgroup(pg_num);
      pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
      get_axes(axes, pointgroup.laue, &pointsym);
      set_transformation_matrix(transform_mat, axes);
    } else {
      pointgroup = ptg_get_pointgroup(0);
    }
  
    return pointgroup;
  }
  
  Pointgroup ptg_get_pointgroup(const int pointgroup_number)
  {
    int i;
    Pointgroup pointgroup;
    PointgroupType pointgroup_type;
  
    pointgroup.number = pointgroup_number;
    pointgroup_type = pointgroup_data[pointgroup_number];
    strcpy(pointgroup.symbol, pointgroup_type.symbol);
    strcpy(pointgroup.schoenflies, pointgroup_type.schoenflies);
    for (i = 0; i < 5; i++) {
      if (pointgroup.symbol[i] == ' ') {pointgroup.symbol[i] = '\0';}
    }
    for (i = 0; i < 3; i++) {
      if (pointgroup.schoenflies[i] == ' ') {pointgroup.schoenflies[i] = '\0';}
    }
    pointgroup.holohedry = pointgroup_type.holohedry;
    pointgroup.laue = pointgroup_type.laue;
  
    debug_print("ptg_get_pointgroup: %s\n", pointgroup_type.symbol);
  
    return pointgroup;
  }
  
  PointSymmetry ptg_get_pointsymmetry(SPGCONST int rotations[][3][3],
                                      const int num_rotations)
  {
    int i, j;
    PointSymmetry pointsym;
  
    pointsym.size = 0;
    for (i = 0; i < num_rotations; i++) {
      for (j = 0; j < pointsym.size; j++) {
        if (mat_check_identity_matrix_i3(rotations[i], pointsym.rot[j])) {
          goto escape;
        }
      }
      mat_copy_matrix_i3(pointsym.rot[pointsym.size], rotations[i]);
      pointsym.size++;
    escape:
      ;
    }
  
    return pointsym;
  }
  
  static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
                                                const int num_rotations)
  {
    PointSymmetry pointsym;
  
    pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
    return get_pointgroup_number(&pointsym);
  }
  
  static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym)
  {
    int i, j, pg_num, counter;
    int table[10];
    PointgroupType pointgroup_type;
  
    debug_print("get_pointgroup_number:");
  
  
    pg_num = 0;
  
    /* Get list of point symmetry operations */
    if (! get_pointgroup_class_table(table, pointsym)) {
      goto end;
    }
  
    for (i = 1; i < 33; i++) {
      counter = 0;
      pointgroup_type = pointgroup_data[i];
      for (j = 0; j < 10; j++) {
        if (pointgroup_type.table[j] == table[j]) {counter++;}
      }
      if (counter == 10) {
        pg_num = i;
        break;
      }
    }
  
   end:
    debug_print(" %d\n", pg_num);
    return pg_num;
  }
  
  static int get_pointgroup_class_table(int table[10],
                                        SPGCONST PointSymmetry * pointsym)
  {
    /* Look-up table */
    /* Operation   -6 -4 -3 -2 -1  1  2  3  4  6 */
    /* Trace     -  2 -1  0  1 -3  3 -1  0  1  2 */
    /* Determinant -1 -1 -1 -1 -1  1  1  1  1  1 */
  
    /* table[0] = -6 axis */
    /* table[1] = -4 axis */
    /* table[2] = -3 axis */
    /* table[3] = -2 axis */
    /* table[4] = -1 axis */
    /* table[5] =  1 axis */
    /* table[6] =  2 axis */
    /* table[7] =  3 axis */
    /* table[8] =  4 axis */
    /* table[9] =  6 axis */
  
    int i, rot_type;
  
    for (i = 0; i < 10; i++) { table[i] = 0; }
    for (i = 0; i < pointsym->size; i++) {
      rot_type = get_rotation_type(pointsym->rot[i]);
      if (rot_type == -1) {
        goto err;
      } else {
        table[rot_type]++;
      }
    }
  
    return 1;
  
   err:
    warning_print("spglib: No point group symbol found ");
    warning_print("(line %d, %s).\n", __LINE__, __FILE__);
    return 0;
  }
  
  static int get_rotation_type(SPGCONST int rot[3][3])
  {
    int rot_type;
  
    if (mat_get_determinant_i3(rot) == -1) {
      switch (mat_get_trace_i3(rot)) {
      case -2: /* -6 */
        rot_type = 0;
        break;
      case -1: /* -4 */
        rot_type = 1;
        break;
      case 0:  /* -3 */
        rot_type = 2;
        break;
      case 1:  /* -2 */
        rot_type = 3;
        break;
      case -3: /* -1 */
        rot_type = 4;
        break;
      default:
        rot_type = -1;
        break;
      }
    } else {
      switch (mat_get_trace_i3(rot)) {
      case 3:  /* 1 */
        rot_type = 5;
        break;
      case -1: /* 2 */
        rot_type = 6;
        break;
      case 0:  /* 3 */
        rot_type = 7;
        break;
      case 1:  /* 4 */
        rot_type = 8;
        break;
      case 2:  /* 6 */
        rot_type = 9;
        break;
      default:
        rot_type = -1;
        break;
      }
    }
  
    return rot_type;
  }
  
  static int get_axes(int axes[3],
                      const Laue laue,
                      SPGCONST PointSymmetry * pointsym)
  {
    switch (laue) {
    case LAUE1:
      axes[0] = 0;
      axes[1] = 1;
      axes[2] = 2;
      break;
    case LAUE2M:
      laue2m(axes, pointsym);
      break;
    case LAUEMMM:
      lauennn(axes, pointsym, 2);
      break;
    case LAUE4M:
      laue_one_axis(axes, pointsym, 4);
      break;
    case LAUE4MMM:
      laue_one_axis(axes, pointsym, 4);
      break;
    case LAUE3:
      laue_one_axis(axes, pointsym, 3);
      break;
    case LAUE3M:
      laue_one_axis(axes, pointsym, 3);
      break;
    case LAUE6M:
      laue_one_axis(axes, pointsym, 3);
      break;
    case LAUE6MMM:
      laue_one_axis(axes, pointsym, 3);
      break;
    case LAUEM3:
      lauennn(axes, pointsym, 2);
      break;
    case LAUEM3M:
      lauennn(axes, pointsym, 4);
      break;
    default:
      break;
    }
  
    return 1;
  }
  
  static int laue2m(int axes[3],
                    SPGCONST PointSymmetry * pointsym)
  {
    int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
    int prop_rot[3][3], t_mat[3][3];
    int ortho_axes[NUM_ROT_AXES];
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search two-fold rotation */
      if (! (mat_get_trace_i3(prop_rot) == -1)) {
        continue;
      }
  
      /* The first axis */
      axes[1] = get_rotation_axis(prop_rot);
      break;
    }
  
    /* The second axis */
    num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 2);
    if (! num_ortho_axis) { goto err; }
  
    min_norm = 8;
    is_found = 0;
    for (i = 0; i < num_ortho_axis; i++) {
      norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
      if (norm < min_norm) {
        min_norm = norm;
        axes[0] = ortho_axes[i];
        is_found = 1;
      }
    }
    if (! is_found) { goto err; }
  
    /* The third axis */
    min_norm = 8;
    is_found = 0;
    for (i = 0; i < num_ortho_axis; i++) {
      norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
      if (norm < min_norm && (! (ortho_axes[i] == axes[0]))) {
        min_norm = norm;
        axes[2] = ortho_axes[i];
        is_found = 1;
      }
    }
    if (! is_found) { goto err; }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[2];
      axes[2] = tmpval;
    }
  
    return 1;
  
   err:
    return 0;
  }
  
  #ifdef SPGDEBUG
  static int lauemmm(int axes[3],
                     SPGCONST PointSymmetry * pointsym)
  {
    int i, count, axis;
    int prop_rot[3][3];
  
  
    for (i = 0; i < 3; i++) { axes[i] = -1; }
    count = 0;
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search two-fold rotation */
      if (! (mat_get_trace_i3(prop_rot) == -1)) {
        continue;
      }
  
      axis = get_rotation_axis(prop_rot);
      if (! ((axis == axes[0]) ||
             (axis == axes[1]) ||
             (axis == axes[2]))) {
        axes[count] = axis;
        count++;
      }
    }
  
    sort_axes(axes);
  
    return 1;
  }
  
  static int laue4m(int axes[3],
                    SPGCONST PointSymmetry * pointsym)
  {
    int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
    int axis_vec[3];
    int prop_rot[3][3], t_mat[3][3];
    int ortho_axes[NUM_ROT_AXES];
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search foud-fold rotation */
      if ( mat_get_trace_i3(prop_rot) == 1) {
        /* The first axis */
        axes[2] = get_rotation_axis(prop_rot);
        break;
      }
    }
  
    /* The second axis */
    num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 4);
    if (! num_ortho_axis) { goto err; }
  
    min_norm = 8;
    is_found = 0;
    for (i = 0; i < num_ortho_axis; i++) {
      norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
      if (norm < min_norm) {
        min_norm = norm;
        axes[0] = ortho_axes[i];
        is_found = 1;
      }
    }
    if (! is_found) { goto err; }
  
    /* The third axis */
    mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
    is_found = 0;
    for (i = 0; i < NUM_ROT_AXES; i++) {
      if (is_exist_axis(axis_vec, i)) {
        is_found = 1;
        axes[1] = i;
        break;
      }
    }
    if (! is_found) { goto err; }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[1];
      axes[1] = tmpval;
    }
  
    return 1;
  
   err:
    return 0;
  }
  
  static int laue4mmm(int axes[3],
                      SPGCONST PointSymmetry * pointsym)
  {
    int i, is_found, tmpval, axis;
    int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
    int axis_vec[3];
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search foud-fold rotation */
      if (mat_get_trace_i3(prop_rot) == 1) {
        /* The first axis */
        axes[2] = get_rotation_axis(prop_rot);
        break;
      }
    }
  
    is_found = 0;
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot2, pointsym->rot[i]);
  
      /* Search two-fold rotation */
      if (! (mat_get_trace_i3(prop_rot2) == -1)) {
        continue;
      }
  
      /* The second axis */
      axis = get_rotation_axis(prop_rot2);
      if (! (axis == axes[2])) {
        axes[0] = axis;
        is_found = 1;
        break;
      }
    }
    if (! is_found) { goto err; }
  
  
    /* The third axis */
    mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
    is_found = 0;
    for (i = 0; i < NUM_ROT_AXES; i++) {
      if (is_exist_axis(axis_vec, i)) {
        is_found = 1;
        axes[1] = i;
        break;
      }
    }
    if (! is_found) { goto err; }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[1];
      axes[1] = tmpval;
    }
  
    return 1;
  
   err:
    return 0;
  }
  
  static int laue3(int axes[3],
                   SPGCONST PointSymmetry * pointsym)
  {
    int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
    int prop_rot[3][3], t_mat[3][3];
    int axis_vec[3];
    int ortho_axes[NUM_ROT_AXES];
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search thee-fold rotation */
      if (mat_get_trace_i3(prop_rot) == 0) {
        /* The first axis */
        axes[2] = get_rotation_axis(prop_rot);
        break;
      }
    }
  
    /* The second axis */
    num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 3);
    if (! num_ortho_axis) { goto err; }
    min_norm = 8;
    is_found = 0;
    for (i = 0; i < num_ortho_axis; i++) {
      norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
      if (norm < min_norm) {
        min_norm = norm;
        axes[0] = ortho_axes[i];
        is_found = 1;
      }
    }
    if (! is_found) { goto err; }
  
    /* The third axis */
    mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
    is_found = 0;
    for (i = 0; i < NUM_ROT_AXES; i++) {
      is_found = is_exist_axis(axis_vec, i);
      if (is_found == 1) {
        axes[1] = i;
        break;
      }
      if (is_found == -1) {
        axes[1] = i + NUM_ROT_AXES;
        break;
      }
    }
    if (! is_found) { goto err; }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[1];
      axes[1] = tmpval;
    }
  
    return 1;
  
   err:
    return 0;
  }
  
  static int laue3m(int axes[3],
                    SPGCONST PointSymmetry * pointsym)
  {
    int i, is_found, tmpval, axis;
    int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
    int axis_vec[3];
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search three-fold rotation */
      if (mat_get_trace_i3(prop_rot) == 0) {
        /* The first axis */
        axes[2] = get_rotation_axis(prop_rot);
        debug_print("laue3m prop_rot\n");
        debug_print_matrix_i3(prop_rot);
        break;
      }
    }
  
    is_found = 0;
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot2, pointsym->rot[i]);
  
      /* Search two-fold rotation */
      if (! (mat_get_trace_i3(prop_rot2) == -1)) {
        continue;
      }
  
      /* The second axis */
      axis = get_rotation_axis(prop_rot2);
      if (! (axis == axes[2])) {
        axes[0] = axis;
        is_found = 1;
        break;
      }
    }
    if (! is_found) { goto err; }
  
    /* The third axis */
    mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
    is_found = 0;
    for (i = 0; i < NUM_ROT_AXES; i++) {
      is_found = is_exist_axis(axis_vec, i);
      if (is_found == 1) {
        axes[1] = i;
        break;
      }
      if (is_found == -1) {
        axes[1] = i + NUM_ROT_AXES;
        break;
      }
    }
    if (! is_found) { goto err; }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[1];
      axes[1] = tmpval;
    }
  
    return 1;
  
   err:
    return 0;
  }
  
  static int lauem3m(int axes[3],
                     SPGCONST PointSymmetry * pointsym)
  {
    int i, count, axis;
    int prop_rot[3][3];
  
    for (i = 0; i < 3; i++) { axes[i] = -1; }
    count = 0;
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search four-fold rotation */
      if (! (mat_get_trace_i3(prop_rot) == 1)) {
        continue;
      }
  
      axis = get_rotation_axis(prop_rot);
      if (! ((axis == axes[0]) ||
             (axis == axes[1]) ||
             (axis == axes[2]))) {
        axes[count] = axis;
        count++;
      }
    }
  
    sort_axes(axes);
  
    return 1;
  }
  #endif
  
  static int laue_one_axis(int axes[3],
                           SPGCONST PointSymmetry * pointsym,
                           const int rot_order)
  {
    int i, j, num_ortho_axis, det, is_found, tmpval;
    int axis_vec[3], tmp_axes[3];
    int prop_rot[3][3], t_mat[3][3];
    int ortho_axes[NUM_ROT_AXES];
  
    debug_print("laue_one_axis with rot_order %d\n", rot_order);
  
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search foud-fold rotation */
      if (rot_order == 4) {
        if (mat_get_trace_i3(prop_rot) == 1) {
          /* The first axis */
          axes[2] = get_rotation_axis(prop_rot);
          break;
        }
      }
  
      /* Search three-fold rotation */
      if (rot_order == 3) {
        if (mat_get_trace_i3(prop_rot) == 0) {
          /* The first axis */
          axes[2] = get_rotation_axis(prop_rot);
          break;
        }
      }
    }
  
    /* Candidates of the second axis */
    num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, rot_order);
    if (! num_ortho_axis) { goto err; }
  
    tmp_axes[1] = -1;
    tmp_axes[2] = axes[2];
    for (i = 0; i < num_ortho_axis; i++) {
      is_found = 0;
      tmp_axes[0] = ortho_axes[i];
      mat_multiply_matrix_vector_i3(axis_vec,
                                    prop_rot,
                                    rot_axes[tmp_axes[0]]);
      for (j = 0; j < num_ortho_axis; j++) {
        is_found = is_exist_axis(axis_vec, ortho_axes[j]);
        if (is_found == 1) {
          tmp_axes[1] = ortho_axes[j];
          break;
        }
        if (is_found == -1) {
          tmp_axes[1] = ortho_axes[j] + NUM_ROT_AXES;
          break;
        }
      }
  
      if (!is_found) { continue; }
  
      set_transformation_matrix(t_mat, tmp_axes);
      det = abs(mat_get_determinant_i3(t_mat));
      if (det < 4) { /* to avoid F-center choice det=4 */
        axes[0] = tmp_axes[0];
        axes[1] = tmp_axes[1];
        goto end;
      }
    }
  
   err: /* axes are not correctly found. */
    warning_print("spglib: Secondary axis is not found.");
    warning_print("(line %d, %s).\n", __LINE__, __FILE__);
    return 0;
  
   end:
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      tmpval = axes[0];
      axes[0] = axes[1];
      axes[1] = tmpval;
    }
  
    debug_print("axes[0] = %d\n", axes[0]);
    debug_print("axes[1] = %d\n", axes[1]);
    debug_print("axes[2] = %d\n", axes[2]);
  
    return 1;
  
  }
  
  static int lauennn(int axes[3],
                     SPGCONST PointSymmetry * pointsym,
                     const int rot_order)
  {
    int i, count, axis;
    int prop_rot[3][3];
  
    for (i = 0; i < 3; i++) {
      axes[i] = -1;
    }
  
    count = 0;
    for (i = 0; i < pointsym->size; i++) {
      get_proper_rotation(prop_rot, pointsym->rot[i]);
  
      /* Search two- or four-fold rotation */
      if ((mat_get_trace_i3(prop_rot) == -1 && rot_order == 2) ||
          (mat_get_trace_i3(prop_rot) == 1 && rot_order == 4)) {
        axis = get_rotation_axis(prop_rot);
        if (! ((axis == axes[0]) ||
               (axis == axes[1]) ||
               (axis == axes[2]))) {
          axes[count] = axis;
          count++;
        }
      }
    }
  
    sort_axes(axes);
  
    return 1;
  }
  
  static int get_rotation_axis(SPGCONST int proper_rot[3][3])
  {
    int i, axis = -1;
    int vec[3];
  
    /* No specific axis for I and -I */
    if (mat_check_identity_matrix_i3(proper_rot, identity)) {
      goto end;
    }
  
    /* Look for eigenvector = rotation axis */
    for (i = 0; i < NUM_ROT_AXES; i++) {
      mat_multiply_matrix_vector_i3(vec, proper_rot, rot_axes[i]);
      if (vec[0] == rot_axes[i][0] &&
          vec[1] == rot_axes[i][1] &&
          vec[2] == rot_axes[i][2]) {
        axis = i;
        break;
      }
    }
  
   end:
  #ifdef SPGDEBUG
    if (axis == -1) {
      printf("rotation axis cound not found.\n");
    }
  #endif
  
    return axis;
  }
  
  static int get_orthogonal_axis(int ortho_axes[],
                                 SPGCONST int proper_rot[3][3],
                                 const int rot_order)
  {
    int i, num_ortho_axis;
    int vec[3];
    int sum_rot[3][3], rot[3][3];
  
    num_ortho_axis = 0;
  
    mat_copy_matrix_i3(sum_rot, identity);
    mat_copy_matrix_i3(rot, identity);
    for (i = 0; i < rot_order - 1; i++) {
      mat_multiply_matrix_i3(rot, proper_rot, rot);
      mat_add_matrix_i3(sum_rot, rot, sum_rot);
    }
  
    for (i = 0; i < NUM_ROT_AXES; i++) {
      mat_multiply_matrix_vector_i3(vec, sum_rot, rot_axes[i]);
      if (vec[0] == 0 &&
          vec[1] == 0 &&
          vec[2] == 0) {
        ortho_axes[num_ortho_axis] = i;
        num_ortho_axis++;
      }
    }
  
    return num_ortho_axis;
  }
  
  static void get_proper_rotation(int prop_rot[3][3],
                                  SPGCONST int rot[3][3])
  {
    if (mat_get_determinant_i3(rot) == -1) {
      mat_multiply_matrix_i3(prop_rot, inversion, rot);
    } else {
      mat_copy_matrix_i3(prop_rot, rot);
    }
  }
  
  static void set_transformation_matrix(int tmat[3][3],
                                        const int axes[3])
  {
    int i, j, s[3];
  
    for (i = 0; i < 3; i++) {
      if (axes[i] < NUM_ROT_AXES) {
        s[i] = 1;
      } else {
        s[i] = -1; /* axes[i] + NUM_ROT_AXES means improper rotation. */
      }
    }
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++) {
        tmat[i][j] = s[j] * rot_axes[axes[j]%NUM_ROT_AXES][i];
      }
    }
  }
  
  static int is_exist_axis(const int axis_vec[3], const int axis_index)
  {
    if ((axis_vec[0] == rot_axes[axis_index][0]) &&
        (axis_vec[1] == rot_axes[axis_index][1]) &&
        (axis_vec[2] == rot_axes[axis_index][2])) { return 1; }
    if ((axis_vec[0] == -rot_axes[axis_index][0]) &&
        (axis_vec[1] == -rot_axes[axis_index][1]) &&
        (axis_vec[2] == -rot_axes[axis_index][2])) { return -1; }
    return 0;
  }
  
  static void sort_axes(int axes[3])
  {
    int axis;
    int t_mat[3][3];
  
    if (axes[1] > axes[2]) {
      axis = axes[1];
      axes[1] = axes[2];
      axes[2] = axis;
    }
  
    if (axes[0] > axes[1]) {
      axis = axes[0];
      axes[0] = axes[1];
      axes[1] = axis;
    }
  
    if (axes[1] > axes[2]) {
      axis = axes[1];
      axes[1] = axes[2];
      axes[2] = axis;
    }
  
    set_transformation_matrix(t_mat, axes);
    if (mat_get_determinant_i3(t_mat) < 0) {
      axis = axes[1];
      axes[1] = axes[2];
      axes[2] = axis;
    }
  }
  