Periodic boundary condition and minimum image in E8 lattice (8-dimensional dense packing)

Dnpluspbc.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 8

/* latticepoint calculates the Dn minimum image */
void latticepoint(double p[DIM], int32_t mirror[DIM])
{
  int rp, gind = 0, summr = 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);
    summr += mirror[rp];
    p[rp] -= mirror[rp];
  }
  if(summr%2 != 0)
  {
    for(rp=0; rp<DIM; ++rp)
    {
      dtmp = fabs(p[rp]);
      if(maxdelta < dtmp)
      {
        maxdelta = dtmp;
        gind = rp;
      }
    }
    // use next nearest neighbor image
    if(p[gind] < 0)
    {
      --mirror[gind];
      ++p[gind];
    }
    else
    {
      ++mirror[gind];
      --p[gind];
    }
  }
}

/**
 * Minimum image decomposition of the vector
 * @f[ p_mathrm{in} = p_\mathrm{out} + \mathrm{mirror} @f]
 * @f$D_n^+@f$ lattice is a coset of two @f$D_n@f$ lattices 
 * with an offset=(1/2) in all coordinates; valid in even dimension
* @f$D_n^{0+}@f$ lattice is a coset of two @f$D_n@f$ lattices 
 * with an offset=(1/2, ..., 1/2, 0); valid in odd dimension
 *
 * @param[in, out] p DIM-dimensional vector of doubles. Get overwrote 
 * @param[out] mirror store the lattice vector. Initialization not necessary
 * 
 * @return int. 0 - mirror is as it is; 2 - mirror should add the offset
 */ 
int latticepoint_dnplus(double p[DIM], int32_t mirror[DIM])
{
  int rp, rst = 0;
  double norm1=0., norm2 = 0.;
  double pb[DIM];
  int32_t mirrorb[DIM];
#if DIM%2==0
  // even dimension
  for(rp=0; rp<DIM; ++rp)
#else
  // odd dimension
  pb[DIM-1] = p[DIM-1];
  for(rp=0; rp<DIM-1; ++rp)
#endif
  {
    pb[rp] = p[rp] - 0.5;
  }
  latticepoint(p, mirror);
  latticepoint(pb, mirrorb);
  for(rp=0; rp<DIM; ++rp)
  {
    norm1 += p[rp]*p[rp];
    norm2 += pb[rp]*pb[rp];
  }
  if(norm2<norm1)
  {
    // copy pb to p
    for(rp=0; rp<DIM; ++rp)
    {
      p[rp] = pb[rp];
      mirror[rp] = mirrorb[rp];
    }
    rst = 2; // shift (1/2)
  }
  return rst;
}