9
$\begingroup$

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 GeneralUtilities`PrintDefinitions, 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:

(Internal`PrecAccur[x] * 2 * Log[10.]) / Pi <= N[x] ^ 2
$\endgroup$
  • $\begingroup$ 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]]. $\endgroup$ – J. M.'s discontentment Aug 3 at 15:05
  • $\begingroup$ Use higher precision. look at N[{#, N[FresnelS[8 + 1*^-28], #] - FresnelS[N[8 + 1*^-28, #]]} & /@ Range[30, 50, 2]] // Grid $\endgroup$ – Bob Hanlon Aug 3 at 15:17
  • 1
    $\begingroup$ Wow, WRI support replied within two hours. They agree it's not behaving correctly. $\endgroup$ – Michael E2 Aug 3 at 21:12
  • $\begingroup$ 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. $\endgroup$ – innaiz Aug 4 at 8:39
6
$\begingroup$

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}

and assembling FresnelS[] from that yields

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
| improve this answer | |
$\endgroup$
  • $\begingroup$ Actually, I just discovered the code for FresnelS[] is available via GeneralUtilities`PrintDefinitions@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.) $\endgroup$ – Michael E2 Aug 3 at 17:38
  • $\begingroup$ It does look to be so when I inspected it myself, and now the fault is not overly surprising. Hopefully it gets resolved soon. $\endgroup$ – J. M.'s discontentment Aug 4 at 1:20
  • 1
    $\begingroup$ 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). $\endgroup$ – Ruslan Aug 4 at 9:43
4
$\begingroup$

The incorrect result is:

FresnelS[N[8+1*^-28, 32]] //InputForm

0.4607524835944079246`3.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.4607524835944079246`3.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.4602142143930144838619886320710524339132605359776113605072`40.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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ +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. $\endgroup$ – Michael E2 Aug 3 at 17:17
-3
$\begingroup$
ListPlot[Table[10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]], {n, 0, 99}]]

ListPlot of the Accuracy of relevance

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}]]

ListPlot from 1 to 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}]]

ListPlot 11 to 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 31 to 45

ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-28, n]]]}, {n, 32, 
   45, 1}]]

ListPlot 32 to 45

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}]]

ListPlot from 44 to 48

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}]]

ListPlot without potentiation 1 to 31

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

GeneralUtilities`PrintDefinitions@FresnelS

...

are much more relevant.

GeneralUtilities`PrintDefinitions@FresnelS

... HoldPattern[FresnelS][ Pattern[TrigExpIntegralDump`x, 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

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

0``54.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}]]

enter image description here

ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-150, n]]]}, {n, 1, 
   46, 1}]]

enter image description here

ListPlot[Table[{n, 10^-Accuracy[FresnelS[N[8 + 1*^-150, n]]]}, {n, 1, 
   46, 1}], PlotRange -> Full]

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}]]

Plot of the difference

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}]]

larger plot

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}]]

Plot

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}]]

limiting behavior

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.