j1 function

double j1(
  1. double x
)

Implementation

double j1(double x) {
  double sign = x < 0 ? -1.0 : 1.0;
  x = x.abs();
  if (x < 8.0) {
    double sum = 0, term = x / 2.0;
    for (int k = 0; k < 60; k++) {
      sum += term;
      term *= -(x * x) / (4.0 * (k + 1) * (k + 2));
      if (term.abs() < 1e-17 * sum.abs() && k > 0) break;
    }
    return sign * sum;
  }
  double theta = x - 3 * M_PI_4;
  return sign * math.sqrt(M_2_PI / x) * math.cos(theta);
}