# The numerics of FresnelS[]

Bug introduced in 12.1 or earlier and persisting through 12.1.1 or later

[CASE:4615361]

Note: A worse problem existed in 12.0 for inputs greater than 8 and of precisions less than 43.66; 12.1 fixed the problem for precisions less than around 32, but the problem for precisions between ~32 and ~43 remains.

The values of N[FresnelS[8 + 1*^-28], 32] and FresnelS[N[8 + 1*^-28, 32]] are surprisingly far apart:

N[FresnelS[8 + 1*^-28], 32] - FresnelS[N[8 + 1*^-28, 32]]
(*  -0.0005  *)


I suppose N[FresnelS[8 + 1*^-28], 32] is the more accurate value, but how can I be sure? What is it's true error?

Here's a broader view:

Block[{$MaxExtraPrecision = 500}, ListLinePlot[ Table[N[FresnelS[x], 32] - FresnelS[N[x, 32]] // RealExponent, {x, Subdivide[0, 15, 15*30]}], PlotRange -> {-36.5, 0.3}, DataRange -> {0, 15}] ]  Is this a bug? Or some inevitable numerical difficulty? How to accurately evaluate FresnelS[x]? Update After @J.M.'s and @Carl's answers, I looked at Trace to see if there were any clue why the accuracy would jump around a precision of 43/44. I discovered the code for the Fresnel family of functions is exposed and can be inspected with GeneralUtilitiesPrintDefinitions, which is unusual for System  functions. I'm not an expert on the Fresnel functions, so sorting it out will take much longer than it's worth to me. I'll happily leave that to the WRI developers. There seems to be a less egregious problem with the machine-precision computation of FresnelS[x] as x increases toward 4, which can be seen in the ramp of the first plot here. The cut-off prec = 43.66 for x = 8 + 1*^-28 in FresnelS[N[x, prec]] actually depends on x and can be found in the code for FresnelF: (InternalPrecAccur[x] * 2 * Log[10.]) / Pi <= N[x] ^ 2  • To repeat the observation I made in the chat room for people who don't follow it: the "bad" evaluation also happens if you try expressing the Fresnel integral in terms of the auxiliary functions, which are (in theory) supposed to be more stable: With[{x = N[8 + 1*^-28, 32]}, 1/2 - FresnelF[x] Cos[Pi x^2/2] - FresnelG[x] Sin[Pi x^2/2]]. – J. M.'s discontentment Aug 3 at 15:05 • Use higher precision. look at N[{#, N[FresnelS[8 + 1*^-28], #] - FresnelS[N[8 + 1*^-28, #]]} & /@ Range[30, 50, 2]] // Grid – Bob Hanlon Aug 3 at 15:17 • Wow, WRI support replied within two hours. They agree it's not behaving correctly. – Michael E2 Aug 3 at 21:12 • Nice. In versions 5.2 and 8.0.4 of Mathematica, N[FresnelS[8 + 1*^-28], 32] - FresnelS[N[8 + 1*^-28, 32]] is simply 0.*^-31. The progress cannot be stopped, however. – innaiz Aug 4 at 8:39 ## 3 Answers I think this is worth reporting to Support. For instance, using formula 7.5.8 from the DLMF: With[{x = N[8 + 1*^-28, 32]}, With[{ζ = Sqrt[π] (1 - I) x/2}, Im[(1 + I)/2 Erf[ζ]]]] 0.46021421439301448386198863207105  and the result is comparable to evaluating N[FresnelS[8 + 1*^-28], 32]. In theory, one is supposed to use the auxiliary functions $$f(z)$$ and $$g(z)$$ for computing Fresnel integrals of moderate or large arguments. However, With[{x = N[8 + 1*^-28, 32]}, {FresnelG[x], FresnelF[x]}] {0.00019781962280286444301613974000765, 0.0392}  With[{x = N[8 + 1*^-28, 32]}, {gg, ff} = {FresnelG[x], FresnelF[x]}; 1/2 - ff Cos[π x^2/2] - gg Sin[π x^2/2]] 0.4608  Contrast this with (cf. formula 7.5.10): With[{x = N[8 + 1*^-28, 32]}, With[{ζ = Sqrt[π] (1 - I) x/2}, ReIm[(1 + I)/2 Exp[ζ^2] Erfc[ζ]]]] {0.000197819622802864443016139740, 0.039785785606985516138011367928}  which works much better: With[{x = N[8 + 1*^-28, 32]}, With[{ζ = Sqrt[π] (1 - I) x/2}, {gg, ff} = ReIm[(1 + I)/2 Exp[ζ^2] Erfc[ζ]]]; 1/2 - ff Cos[Pi x^2/2] - gg Sin[Pi x^2/2]] 0.460214214393014483861988632071  • Actually, I just discovered the code for FresnelS[] is available via GeneralUtilitiesPrintDefinitions@FresnelS (or Trace). It does just what you show. You can also inspect FresnelF. (I think there's a problem with FresnelF and maybe G.) – Michael E2 Aug 3 at 17:38 • It does look to be so when I inspected it myself, and now the fault is not overly surprising. Hopefully it gets resolved soon. – J. M.'s discontentment Aug 4 at 1:20 • I wonder whether special functions are even something of a priority to WRI, given that the broken Mathieu functions remain broken for years (even after reporting). – Ruslan Aug 4 at 9:43 The incorrect result is: FresnelS[N[8+1*^-28, 32]] //InputForm  0.46075248359440792463.970167826243401 Note that the precision of the output is 3.97, indicating that the 4th digit may not be accurate, which is exactly what you observe. If you increase the precision: FresnelS[N[8+1*^-28, 43]] //InputForm  0.46075248359440792463.970167826243401 You still only get approximately 4 digits of precision. A further increase to 44 does produce a better answer: FresnelS[N[8+1*^-28, 44]] //InputForm  0.460214214393014483861988632071052433913260535977611360507240.31143196057397 I think having the same result for 32 and 43 digits of precision looks suspicious, and the fact that approximately an additional 16 digits of precision (beyond the 28) are needed to get a correct result suggests that there may be some incorrect machine number approximation being used under the hood. I would also suggest reporting this issue to support. • +1 & thanks. I will point out that the uncertainty 10^-Accuracy[FresnelS[N[8 + 1*^-28, 32]]] is a bit less than the error, so I think the precision 3.97 is higher than it should be. The constant result for precision between 32 and 43 is suspicious. I'll report it. – Michael E2 Aug 3 at 17:17 ListPlot[Table[10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]], {n, 0, 99}]]  This alters the cognition of what is really going on hopefully. This plot causes some more sensation: ListPlot[Table[10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]], {n, 0, 10}]]  How can this happen? It is the very same function on a subinterval. ListPlot[Table[10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]], {n, 11, 99}]]  More as expected the Accuracy falls monotonically but slow. ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]]}, {n, 31, 45, 1}]]  ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]]}, {n, 32, 45, 1}]]  10^-Accuracy[FresnelS[N[8 + 1*^-28, 32]]]  0.0000493514 ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]]}, {n, 44, 48, 1}]]  This is rapid decreasing and not some magnitudes smaller but 10^-Accuracy[FresnelS[N[8 + 1*^-28, 32]]]/ 10^-Accuracy[FresnelS[N[8 + 1*^-28, 44]]] 2.19671*10^36  36 magnitudes smaller. ListPlot[Table[{n, Accuracy[FresnelS[N[8 + 1*^-28, n]]]}, {n, 1, 31, 1}]]  FresnelS[] is on the summary of new features of version 12. So this can not be a bug. It is intended. The annotation is "(updated) — substantially expanded examples for 250 special functions". On the page linked from there, Guide: Mathematical Functions they mention "The Wolfram Language has the most extensive collection of mathematical functions ever assembled. Often relying on original results and algorithms developed at Wolfram Research over the past two decades, each function supports a full range of symbolic operations, as well as efficient numerical evaluation to arbitrary precision, for all complex values of parameters." So this behavior gives a look into the efficient numerical evaluation to arbitrary precision of the Wolfram Language. MathematicalFunctionData["FresnelS", "RelatedFunctionRepresentations"] {Function[{\[FormalZ]}, Inactivate[ FresnelS[\[FormalZ]] == 1/4 (1 + I) (Erf[1/2 (1 + I) Sqrt[\[Pi]] \[FormalZ]] - I Erf[1/2 (1 - I) Sqrt[\[Pi]] \[FormalZ]])]], Function[{\[FormalZ]}, Inactivate[ FresnelS[\[FormalZ]] == 1/4 (I - 1) (Erfi[1/2 (1 - I) Sqrt[\[Pi]] \[FormalZ]] + I Erfi[1/2 (1 + I) Sqrt[\[Pi]] \[FormalZ]])]], Function[{\[FormalZ]}, Inactivate[ FresnelS[\[FormalZ]] == 1/2 + 1/4 (1 + I) (-Erfc[1/2 (1 + I) Sqrt[\[Pi]] \[FormalZ]] + I Erfc[1/2 (1 - I) Sqrt[\[Pi]] \[FormalZ]])]], Function[{\[FormalZ]}, Inactivate[ FresnelS[\[FormalZ]] == ((1 - I) (DawsonF[1/2 (I - 1) Sqrt[\[Pi]] \[FormalZ]] - I E^(I \[Pi] \[FormalZ]^2) DawsonF[ 1/2 (1 + I) Sqrt[\[Pi]] \[FormalZ]]))/((2 Sqrt[\[Pi]]) E^( 1/2 I \[Pi] \[FormalZ]^2))]], Function[{\[FormalZ]}, Inactivate[ FresnelS[\[FormalZ]] == FresnelF[\[FormalZ]] (-Cos[(\[Pi] \[FormalZ]^2)/2]) - FresnelG[\[FormalZ]] Sin[(\[Pi] \[FormalZ]^2)/2] + 1/2]]}  So the behavior is brand new and at the kernel of the innovations. The page MathematicalFunctionData states all definitions are now direct Digital Library of Mathematics Functions standard. MathematicalFunctionData["FresnelS", "ArgumentPattern"] Inactive[FresnelS][_] MathematicalFunctionData["FresnelS", "WolframFunctionsSiteLink"] Hyperlink["http://functions.wolfram.com/GammaBetaErf/FresnelS/", \ "http://functions.wolfram.com/GammaBetaErf/FresnelS/"]  This shows another hierarchical depends used internally to calculate the values the accuracy is derived from. The overwhelming power of mathematical knowledge representation leads to the Mathematica input: MathematicalFunctionData["FresnelS", "PropertyAssociation"] ...  For numerical evaluation the definitions in GeneralUtilitiesPrintDefinitions@FresnelS ...  are much more relevant. GeneralUtilitiesPrintDefinitions@FresnelS ... HoldPattern[FresnelS][ Pattern[TrigExpIntegralDumpx, Blank[Real]]] := Module[{z = Pi * x ^ 2, ax}, If[Less[z, 8.], Times[(z * x) / 6, HypergeometricPFQ[{3 / 4}, {3 / 2, 7 / 4}, -(z / 4) ^ 2 ] ], DivideBy[z, 2]; ax = Abs @ x; Times[Sign @ x, (1 / 2) - ((FresnelG[ax] * Sin[z]) + FresnelF[ax] * Cos[z]) ] ] ]; ... This has a conditional at z=8.. At this point on the Reals representation to calculate the numerical values is changed. That holds for the jump in Accuracy. fslist = {FresnelS[N[8 - 1*^-43 + 8 1*^-44, 100]] // InputForm, FresnelS[N[8 - 1*^-43 + 9 1*^-44, 100]] // InputForm, FresnelS[N[8 - 1*^-43 + 1*^-44, 100]] // InputForm, FresnelS[N[8 - 1*^-43 + 11 1*^-44, 100]] // InputForm, FresnelS[N[8 - 1*^-43 + 12 1*^-44, 100]] // InputForm} // Dataset  There is a jump in accuracy and in the value at 8. To get this Dataset it is necessary to use the 8 and the notation 1*^-43 otherwise the represetion is not that exact or change the corresponding values in the notebooks settings. (FresnelS[N[8 - 1*^-43 + 1*^-44, 100]] - FresnelS[N[8 - 1*^-43 + 11 1*^-44, 100]]) // InputForm 054.488171823879036  A very accurate zero. But not machine presicion or arbitrary precision. There is evidence that the considerations of the designer were reduce the time for calculations first then presicion at 8 the situation changes and the higher precision approximation for lower values is replaced by the representation with higher precision. Since all other have less precision this suffices to be market leader again disregading Fortran implementation for arbitrary presicion. Since these values can not be exact to arbitrary presision but just to a very high precision this choice was there and is acceptable. I get a bigger problem with this ListPlot[Table[{n, Accuracy[FresnelS[N[8 + 1*^-150, n]]]}, {n, 1, 46, 1}]]  ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-150, n]]]}, {n, 1, 46, 1}]]  ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-150, n]]]}, {n, 1, 46, 1}], PlotRange -> Full]  So the accuracy get smaller the higher the value for Accurcay is selected. It has a clear maximum at the Accuracy n=4 at 5 10^41 and crosses zero between n=45 and n=46. Expected is an Accuracy that is negative and gets more negative. This is a behaviour forced by the function representation to the Gaussian error function with pure complex arguments. Per definitionem this is so. So to work with FresnelS it is advisable to set the Accuracy parameter bigger than 46 to calculate with sensible values. So the plot under consideration looks different: Block[{$MaxExtraPrecision = 500},
ListLinePlot[
Table[N[FresnelS[x], 50] - FresnelS[N[x, 50]] // RealExponent, {x,
Subdivide[0, 15, 15*30]}], PlotRange -> {-60, 0.3},
DataRange -> {0, 15}]]


This shows there is a strategy necessary to calculate the values the function reaches -50 almost regularly. It is never worse than -3.25.

Block[{$MaxExtraPrecision = 500}, ListLinePlot[ Table[N[FresnelS[x], 50] - FresnelS[N[x, 50]] // RealExponent, {x, Subdivide[0, 15, 15*30]}], PlotRange -> {-56, -3.25}, DataRange -> {0, 15}]]  This about -3 is the difference that almost reached with the -0.0005 in the question. The behavior is even more brave for Block[{$MaxExtraPrecision = 500},
ListLinePlot[
Table[N[FresnelS[x], 100] - FresnelS[N[x, 100]] // RealExponent, {x,
Subdivide[0, 15, 15*30]}], PlotRange -> {-108, -3.5},
DataRange -> {0, 15}]]


So the advice can go further to set the Accuracy parameter as high as possible to avoid the divergence with upper limit and set the algorithm in operation as seldom as possible.

The x=8 remains for all Accuracy parameters this way but the extent of the first value at which the curve comes close to about -3.5 can be extended. The relation seems to be something like square root between 50 and 100. The minimium is always taken for very small x and reaches if the correction algorithm is set into operation to calculate the values.

Block[{\$MaxExtraPrecision = 500},
ListLinePlot[
Table[N[FresnelS[x], 500] - FresnelS[N[x, 500]] // RealExponent, {x,
Subdivide[0, 15, 15*30]}](*,PlotRange\[Rule]{-108,-3.5}*),
DataRange -> {0, 15}]]


I wish at this point I am able to look deeper into the details, but higher value in N is not better is required for high demands. This even gets the values at the peak very low.