Module fresnel 

The Fresnel integrals C(x) = int_0^x cos(pi t^2/2) dt and S(x) = int_0^x sin(pi t^2/2) dt.

Two regions

Below a crossover (2.5265 at both precisions) the integrands are summed directly: C(x)/x and S(x)/x^3 are both smooth functions of w = x^4, fitted as Chebyshev series and summed by Clenshaw. Above it the standard auxiliary form

C = 1/2 + f sin t - g cos t,   S = 1/2 - f cos t - g sin t,   t = pi x^2 / 2

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

Why Chebyshev below and Horner above

Measurement, not symmetry. The small-argument fit in the monomial basis has an error amplification (sum |c_k| |w|^k / |f|) of 3482 at this crossover, and 1.26e5 if the crossover moves to 3. The same fit in the Chebyshev basis summed by Clenshaw sits at 5.8 and 7.0. Monomial Horner would cap the f64 kernel at about 4e-13. Clenshaw costs two operations per term against Horner’s one and buys three orders of magnitude, which is also what lets the crossover sit far enough out for the auxiliaries to be well conditioned. Their own amplification is 1.02 there, so they keep Horner.

The phase

t = pi x^2 / 2 computed as x*x*0.5 is worthless long before the function is: measured against a 45-digit oracle, sin(pi*(x*x*0.5)) in binary64 is 5.2e-13 off at x = 123, 9.8e-11 at 1234, 5.3e-6 at 98765, and returns the wrong sign by x ~ 1e9. Since C and S are 1/2 plus a term of size 1/(pi x), that error lands directly on the result.

phase_half_x2 fixes it in about ten operations, and the fix is exact: x*x splits as p + e with e always representable, halving is exact, and each half reduces mod 2 exactly by Sterbenz. Both halves must be reduced before being added. |e/2| reaches ulp(x^2)/4, which is 32 at x = 1e9, and adding that to an already-reduced p/2 rounds the latter’s low bits straight off. Measured 2.80 ulp (C) and 2.64 (S) in f64 over x from 1e-4 to 1e15, and 2.14 / 3.40 in f32 out to 1e7. The naive phase alone is worth thousands of ulp there.

Above x = 1.147e16 (f64) / 2.136e7 (f32) the correction has fallen under half an ulp of 1/2 and both functions are exactly 1/2. Below that x^2 cannot overflow, so the phase needs no range guard.

Functions

fresnel_with

(S(x), C(x)), in SciPy’s order.

phase_half_x2

x^2/2 mod 2, the argument for sincos_pi, to full precision for every x whose square is finite.