phd.c

double covariance(double *x, int len_x, int k)
{
  int i;
  double dot = 0;
  for (i = 0; i < len_x - k; ++i){
    dot += x[i] * x[i+k];
  }
  return dot / len_x;
}

double phd(double *x, int len_x)
{
  double r1 = covariance(x, len_x, 1);
  double r2 = covariance(x, len_x, 2);
  double a = (r2 + sqrt(r2 * r2 + 8 * r1 * r1)) / 4 / r1;
  return arccos(a);
}