Periodic boundary condition and minimum image in D4 lattice (4-dimensional dense packing)

Dnpbc.hpp

/* Periodic boundary condition and minimum image in Dn lattice (n-dimensional dense packing, n>=4) 
 * 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>

// Dimensionality DIM >= 4
#define DIM 4

/**
 * Minimum image decomposition of a vector under Dn periodic boundary condition
 * @f[ p_mathrm{in} = p_\mathrm{out} + \mathrm{mirror} @f]
 * For @f$D_n@f$ lattices, @f$\sum_i \mathrm{mirror}_i \equiv 0 \mod 2@f$
 *
 * @param[in, out] p DIM-dimensional vector of doubles. Get overwrote 
 * @param[out] mirror store the lattice vector. Initialization not necessary
 * 
 */ 
void 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];
    }
  }
}