[NAThompson]

fresnel_c, fresnel_s and fresnel_cis (also in Python), within an ulp in float, double and multiprecision:
a double-word Maclaurin series, the even Laplace continued fraction for f and g, and an exact phase.
Also names the temporaries of two_sum_error, which expression templates reassociated.

C(x) = ∫₀ˣ cos(πt²/2) dt and S(x) = ∫₀ˣ sin(πt²/2) dt, the convention of DLMF 7.2, A&S, Mathematica, mpmath, SciPy and Arb (normalized = 1). fresnel_cis(x) returns C(x) + i S(x) = ∫₀ˣ cis(πt²/2) dt for the cost of one of them, as std::complex<T> for float, double and long double, and as the matching Boost.Multiprecision complex type otherwise (via tools::detail::complex_from_real_t). In Python they are boost.math.fresnel_c, fresnel_s and fresnel_cis (complex128), elementwise over arrays like the other special functions.

With θ = πx²/2 and p bits:

• θ ≤ p/10: the Maclaurin series of C + iS, Horner in θ² in double-word arithmetic.
• θ > p/10: C = ½ + f sin θ − g cos θ, S = ½ − f cos θ − g sin θ, with g + i f from the even part of Laplace's continued fraction for erfc, evaluated backwards, last step in double-word.
• x > ε^(−1/2): f = 1/(πx), g = 1/(π²x³).
• The phase: x² = p + e exactly, p/2 reduced mod 2 exactly, e/2 applied by angle addition.

Relative error in ε against mpmath at 80 digits, max / mean:

type	series	continued fraction	asymptotic
float	0.50 / 0.19	0.50 / 0.18	0.08
float, promote_float<false>	0.50 / 0.19	0.75 / 0.19	0.50 / 0.14
double	0.50 / 0.18	0.69 / 0.19	0.50 / 0.19
cpp_bin_float_50	0.49 / 0.18	0.64 / 0.19

mpfr_float is correctly rounded at 50, 100 and 200 digits, also after changing the precision at runtime (precision-dependent tables are cached only for compile-time precision).

google/benchmark, double, -O2, Apple M-series:

path	time
series, x ∈ [0, 1.8]	159 ns
continued fraction, x ∈ [1.9, 3]	277 ns
continued fraction, x ∈ [3, 100]	82 ns
asymptotic	37 ns
fresnel_cis, series	160 ns

two_sum_error in compensated_arithmetic.hpp: with expression templates on (e.g. mpfr_float), (a - (s - bb)) + (b - bb) was reassociated, and for |a| < |b| the error came out as 0 instead of the rounding error of a + b (Fresnel at 100 digits was off by ~1e7 ε). Each difference now has a named temporary; fast_two_sum and the Veltkamp split get the same treatment. Struve, exp_sinh, the Gauss quadratures, Leja, 0F1, 2F0, Laguerre, Kolmogorov–Smirnov, Jacobi theta, beta and ibeta tests pass.

🤖 Generated with Claude Code