Function y0_near_zero 

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

Generated by SageMath: Evaluates: Y0(x) = 2/pi*(euler_gamma + log(x/2))J0(x) - sum((-1)^m(x/2)^(2m)/(m!)^2sum(1+1/2 + … 1/m)) expressed as: Y0(x)=log(x)*W0(x) - Z0(x)

from sage.all import *

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

# Define J0(x) Taylor expansion at x = 0
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)])

J0_series = j_series(0, x)

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

W0 = (d2/pi) * J0_series
Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)

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