Function y1_near_zero 

fn y1_near_zero(x: f64, w_log: DoubleDouble) -> f64

Generated by SageMath: Evaluates: y2 = -J1(x)log(x) + 1/x * (1 - sum((-1)^m(H(m)+H(m-1))/(2^mm!(m-1)!)x^(2m)) Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x)) expressed as: Y1(x)=log(x)W1(x) - Z1(x) - 2/(pix)

from sage.all import *

R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
x = R.gen()
N = 16  # Number of terms (adjust as needed)
gamma = RealField(300)(euler_gamma)
d2 = RealField(300)(2)
pi = RealField(300).pi()
log2 = RealField(300)(2).log()

def j_series(n, x):
    return sum([(-1)**m * (x/2)**(ZZ(n) + ZZ(2)*ZZ(m)) / (ZZ(m).factorial() * (ZZ(m) + ZZ(n)).factorial()) for m in range(N)])

J1_series = j_series(1, x)

def harmony(m):
    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))

def z_series(x):
    return sum([(-1)**m * (x)**(ZZ(2)*ZZ(m)) / (ZZ(2)**(2*m) * ZZ(m).factorial() * (ZZ(m) - ZZ(1)).factorial()) * (harmony(m) + harmony(m - 1)) for m in range(1, N)])

W1 = d2/pi * J1_series
Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))

def y1_full(x):
    return d2/pi * (J1_series(x) * x.log() - 1/x * ( 1 - z_series(x)) + (gamma - log2) * J1_series(x))

# see the series
print(W0)
print(Z0)

See ./notes/bessel_y1_taylor.ipynb for generation