Module sici 

The trigonometric integrals Si(x) = int_0^x sin(t)/t dt and Ci(x) = gamma + ln x + int_0^x (cos t - 1)/t dt.

Two regions

Below a crossover (12 in f64, 6 in f32) Si(x)/x and Cin(x)/x^2 are smooth functions of v = x^2, fitted as Chebyshev series and summed by Clenshaw, where Cin = gamma + ln x - Ci is the entire part, the piece that is not the logarithmic singularity. Above it, the auxiliary form

Si = pi/2 - f cos x - g sin x,   Ci = f sin x - g cos x

with f = P(v)/x and g = Q(v)/x^2 for v = 1/x^2, both Horner polynomials tending to 1.

Why the crossover is so far out

Because these auxiliaries are harder than the Fresnel ones, and for a structural reason worth recording. f(x) = int_0^inf e^{-xt}/(1+t^2) dt (verified to 12 digits), so P(v) = x f is a Stieltjes function: its asymptotic series sum (-1)^k (2k)! v^k diverges, and the branch cut reaches v = 0. Polynomial convergence at that endpoint is therefore sub-geometric. It shows. Measured degree for f64, at contribution-weighted targets:

range	deg P	amplification	deg Q
x >= 6	33	3.6e8	31
x >= 8	25	455	24
x >= 10	20	1.19	19
x >= 12	17	1.03	16

One fit at x >= 12 beats the multi-range ladders that were also measured ([8,20] plus [20,inf) is 28 terms and a select, against 17), and Pade of the divergent series is no better. [10/10] reaches only 8.8e-8 at x >= 8.

P is fitted at plain relative accuracy because |Ci| ~ f, so f’s error is the result’s error. Q is relaxed by a factor x, contributing at 1/x^2 against a 1/x result. Re-fitting Q at plain relative accuracy adds degrees and buys nothing.

Accuracy

Measured against mpmath at 45 digits over x from 1e-4 to 1e15 (f64) and 1e7 (f32): Si 2.03 ulp f64 / 1.34 f32, Ci 1.42 / 1.99 relative to its envelope.

Two contract points belong in the caller’s head:

• Ci has zeros, the first near x = 0.6165, and nothing is relatively accurate at one. The grading above is against |gamma + ln x| + |Cin| below the crossover and 1/x above it, which is what the arithmetic can actually deliver.
• Large-x accuracy inherits sin_cos’s argument reduction. For Ci the oscillation is the value, so a phase error is a relative error. Full reduction is a Best-tier property in this library, and below that Ci’s accuracy at large x degrades with it. Si is insulated: it tends to pi/2 and the oscillation is a correction of size 1/x.

Si is pi/2 to within half an ulp above x = 1.147e16 (f64) / 2.136e7 (f32). Ci has no such cutoff: it decays like 1/x and stays representable for every finite argument.

Functions

sici_with

(Si(x), Ci(x)).