Measuring Entanglement - Code

Real tangle( unsigned int k) const {

  Complex prod = 0; Real A = 0, B = 0, r = 0;

  unsigned long i = 0, j, k1 = 1LU << k; // k1 = 0...010...0, k'th bit = 1

  while( i < N) { // i will have k'th bit = 0, j will have k'th bit = 1

    j = i | k1; prod += Conj(a[i])*a[j];

    A += abs2(a[i]); B += abs2(a[j]);

    ++i; i += (i & k1); // skip i values with k'th bit = 1
  }

  Real t = 2*eps*eps*N; // threshold for zero norm squared: 2*sum{|e+i*e|**2} = 2*(2*e*e*(N/2))

  if( A > t && B > t) r = 1 - abs2(prod)/(A*B);

  return r;
}