/* Periodic boundary condition and minimum image in D4 lattice (4-dimensional dense packing)
 * Created by Yi Hu on 6/13/19.
 * Refer to: J. Convay and N. Sloane, IEEE Trans Inform Theory 28, 227 (1982).
 */

#include <cstdint>
#include <cmath>

#define DIM 4

/* input value:
   p - vector of coordinates. Modify to its minimum image
   mirror - store the lattice point
   return value: 0 (suggest mirror is zeros and p is unchanged)
                 or 1 (mirror is not all zeros and p has been changed)
 */
int latticepoint(double p[DIM], int32_t mirror[DIM])
{
  int rp, gind = 0, summr = 0, is_mird = 0;
  double maxdelta = 0, dtmp;
  for(rp=0; rp<DIM; ++rp)
  {
    // fast round eqv mirror[rp] = round(p[rp]);
    dtmp = p[rp] + 6755399441055744.0;
    mirror[rp] = reinterpret_cast<int32_t&>(dtmp);
    p[rp] -= mirror[rp];
    dtmp = fabs(p[rp]);
    if(maxdelta < dtmp)
    {
      maxdelta = dtmp;
      gind = rp;
    }
    summr += mirror[rp];
  }
  if(summr%2 != 0)
  {
    // use next nearest neighbor image
    if(p[gind] < 0)
    {
      --mirror[gind];
      ++p[gind];
    }
    else
    {
      ++mirror[gind];
      --p[gind];
    }
  }
  for(rp=0; rp<DIM; ++rp)
  {
    if(mirror[rp] != 0)
    {
      is_mird = 1;
      break;
    }
  }
  return is_mird;
}