y0 function

double y0(
  1. double x
)

Implementation

double y0(double x) {
  if (x <= 0) return x == 0 ? -double.infinity : double.nan;
  if (x < 8.0) {
    const euler = 0.5772156649015329;
    double j0val = j0(x);
    double lnx2 = math.log(x / 2);

    double psum = 0;
    double pk = 1.0;
    double hval = 0;
    for (int k = 1; k < 60; k++) {
      hval += 1.0 / k;
      pk *= -(x * x) / (4.0 * k * k);
      psum += -pk * hval;
      if (pk.abs() * hval < 1e-17 * psum.abs() && k > 5) break;
    }
    return (2 / M_PI) * ((euler + lnx2) * j0val + psum);
  }
  double theta = x - M_PI_4;
  return math.sqrt(M_2_PI / x) * math.sin(theta);
}