How do I numerically evaluate and plot the Fabius function?

Question

The Fabius function is a well-known example in analysis of a non-analytic function that is infinitely differentiable. I want to be able to numerically evaluate the function for any real argument, as well as to plot it, just like the picture in Wikipedia:

I had found this implementation by a Mr. Reshetnikov, but it seems to be only for fractions whose denominators are a power of two (a.k.a. dyadic rationals). Can that implementation be extended to compute, say the value at a non-dyadic rational like $2/3$, or a not exact number like $0.775$? Thanks for any ideas.

Answer 1

Given $d\in\mathbb{N}_0$, the Taylor series about $i/2^d$ is a polynomial of degree at most $d$ for all $i\in\mathbb{Z}$.
Let $S_d$ be the set of such Taylor series.

There exist unique polynomials $\text{pol}_0,\text{pol}_1,\ \dots\ ,\text{pol}_d$ of degree $0,1,\ \dots\ ,d$ and a function $c:\mathbb{N}\times\mathbb{R}\mapsto\{-1,0,1\}$, such that for all $x\in\mathbb{R}$, the following sum coincides with those Taylor series in $S_d$ whose point of expansion differs the least from $x$.

$$\sum_{i\ =\ 0}^d c(i, x)\ \text{pol}_i(\Lambda(2-2^{i+1}(x\bmod 2^{1-i})))$$

where $\Lambda$ is HeavisideLambda. The choice of $c$ is unique, except when $x=i/2^{d+1}$, where $i$ is an odd integer. In that case, there are 2 choices because the nearest Taylor series is not unique.

In my approximation, I use that sum of polynomials for computing values on $[0,2)$, and self-similarity for the extension to $\mathbb{R}$ and the derivatives.

Ancient code, see update below:

Young already linked the Mathematics Stack Exchange question with my answer to this question, but I made some changes to the code. The first part which computes the values of $F(2^{-i})$ for $i=1,3,5,\dots$. I have got that to run much faster by avoiding additions of awkward rationals with huge GCDs. With help from OEIS, I found a "Mathematica-closed form" for the needed GCDs, which is much more efficient than the implicit GCD calls from addition of rationals.

Remove[d, recTable, fabs, pol, fabius, y, n, x]

d = 399;   (*polynomial degree of pieces; odd values make fabius[x] continuous*)

recTable = RecurrenceTable;
SetAttributes[{MultiplicativeOrder, pol}, Listable]

fabs = With[{F = Floor[d/2] + 1}, With[{gcd = 2^(F - DigitCount[F, 2, 1]) Times @@ (
         #^Quotient[F, (# - 1)/2]) &[Pick[#, # - 2 MultiplicativeOrder[4, #], 1] &[
           Prime[Range[3, Max[Floor[385/213 #], 4]]]][[;; #]] &[Count[MultiplicativeOrder[
             4, 2 # + 3] - #, 1] &[Range[F - 1]]]]},
  Reap[Fold[Append[#, Sow[Total[#]]] &[# #2] &,
   {Sow[gcd]}, With[{js = Map[# (2 # - 1) &, Range[-F, -2]],
   T = recTable[{y[n + 1] == 4 y[n] + 9, y[1] == 0}, y, {n, F}]},
   Table[Append[T[[i]]/js[[1 - i ;;]], 1], {i, F}]]]][[2, 1]]/FoldList[Times, 2 gcd,
   recTable[{y[n + 2] == 80 y[n + 1] - 1024 y[n], y[1] == 144, y[2] == 11520}, y, {n, F}]]]];

I also changed the form of polynomials, by putting GCDs outside parentheses to avoid repetitive, extensive GCD calculations when calling fabius[x] for different values of x.

Clear[pol]
Evaluate[Table[pol[r, x_], {r, d + 2}]] = With[{F = Floor[d/2]},
 With[{gcd = Prepend[Riffle[#, #2 #], 1][[;; d + 2]]}, 1/gcd
  MapIndexed[If[OddQ[First[#2]], x^-1, x^-2] Fold[(# + #2) x^2 &, 0, #] &,
   With[{p = Riffle[#2/2, #2, {2, d + 1, 2}], S = Join[Range[d, 1, -1], {1}],
    rat = Prepend[Riffle[#2, Ratios[#]/Most[#2]], 1][[;; d + 1]]},
     FoldList[rat[[#2]]/S[[-#2 ;; ;; 2]] If[OddQ[#2], Append[#/p[[#2]],
      2 fabs[[(#2 + 1)/2]] gcd[[#2]]], #/p[[#2]]] &, {}, Range[d + 1]]]]] &[
       Join[{1, 144}, Map[1/GCD @@ (2 fabs[[{3, # - 2, #}]]) &, Range[3, F + 1]]],
       NestList[4 # &, 2, Max[1, F]]]];

The changes to the definition part is hardly significant.

(T0 = Power[-1, #[[;; Ceiling[d/2 + 1]]]]; T1 = Mod[Most[#], 2];
 T2 = Power[2, #]; T3 = Power[2, Rest[#]]) &[Range[0, d + 1]];
Derivative[n_][fabius] = 2^(n (n + 1)/2) fabius[2^n #] &;
With[{prec = N[d/Log2[10] + 15], d2 = d + 2}, fabius[x_?NumericQ] :=
 With[{t = SetPrecision[#, Max[Precision[#], prec]] &[Mod[x, 2]]},
  With[{refs = Position[IntegerDigits[BitXor[BitShiftRight[#, 1], #] &[
   FromDigits[First[RealDigits[t, 2, d2, 0]], 2]], 2, d2], 1][[All, 1]]},
   (-1)^DigitCount[Floor[x, 2], 2, 1] If[
   refs =!= {}, Dot[T0[[Accumulate[Prepend[T1[[Most[refs]]], 1]]]],
  pol[refs, FoldList[(2 - #) #2 &, T2[[First[refs]]] t, T3[[Differences[refs]]]] - 1]], 0]]]]

There you go:

Plot[fabius[x], {x, 0, 24}, AspectRatio -> 1/12, ImageSize -> 800]

If X is uniformly distributed on the unit interval, the expected absolute error and the maximal absolute error of fabius[X] are both rationals. If you use p = ∞ in the following code, you obtain their exact values. All local maxima of the error equal the maximal error, so the approximation is uniform in an absolute sense.

errors = With[{p = 10}, With[{sL = 2 Last[fabs] + Dot[
          Riffle[Table[pol[i + 1, N[1, p + (93 + 5 d)/28]],
           {i, 1, d + 1, 2}], Most[fabs]], Reverse[
            With[{q = Floor[d, 2] + 2}, FoldList[Times, 2^(-q),
             Reverse[recTable[{y[n + 1] == (n - q - 1) y[n]/2/
              (n - q), y[1] == -1/2/q}, y, {n, q - 1}]]]]]]},
               Block[{$MaxExtraPrecision = ∞}, N[#, Precision[#]] &[
                If[OddQ[d], {2^(d + 1) Last[fabs], sL}, {2^(d + 1) sL, fabs[[-2]]}]]]]]

Update:

I'm still using the sum at the top for the approximation, but some functionality is added, and now non-rational input is handled efficiently. Hence I wanted to provide this update:

The form of the polynomials is changed to:
$\ $132809 + 3825 x^2 (2332 + 63 x^2 (266 + 15 x^2 (28 + 9 x^2)))

instead of the HornerForm:
$\ $132809 + x^2 (8919900 + x^2 (64099350 + x^2 (101209500 + 32531625 x^2)))

because the asymptotic memory usage for the latter was very bad. Also, only coefficients for a finitely many (memorizeDegree) polynomial are memoized. Higher degree (maxDegree) polynomials are used, but only one/two of them is in memory at a given time.

I added an argument such that fabius[x, p] returns precision p. Each polynomial is computed to the needed precision only, unless x is rational and ByteCount[x] < 3000. In that case it's simply faster to use rational arithmetic. Such rational input evaluates much faster: E.g. if maxDegree = 1119; memorizeDegree = 449; then Timing[fabius[4/3, 188888];] and Timing[fabius[E, 64888];] are about the same.

fabius[x, ∞] returns the exact result, if possible (i.e., if x is dyadic with sufficiently small denominator). Otherwise, there's a warning.

The initialization part:

Remove[maxDegree, memorizeDegree, track, extraPrecision, fabius,
       polPrecs, polNext, polMove, polCoefs, polEval, logSizes, maxAbsErr]

maxDegree = 999; (* Nonnegativ integer *)
memorizeDegree = 449; (* Integer greater or equal 2; At 1139 my RAM is all used;
                         RAM use increases quadraticly *)
extraPrecision = 4; (* Nonnegativ real *)

polNext[k_, n_] := Module[{L = k},
  L[[1, All, 1]] /= Range[n - 2, 1, -2];
  L[[1]] = Reap[L[[2]] *= 2^n/Fold[Function[G, Sow[#/G]; G][
                 GCD @@ #] &[#2 {1, #}] &, 1/n, First[L]]][[2, 1]];
  If[EvenQ[n], AppendTo[L[[1]], {Fold[First[#2] + Last[#2] # &, 1/(n + 1),
       Transpose[{L[[1, All, 1]]/Range[n - 1, 3, -2], L[[1, All, 2]]}]]/(2^n - 1), 1}];
       (L[[1, -1]] *= #; L[[2]] *= #;) &[Denominator[L[[1, -1, 1]]]]]; L]


Dynamic[ToString[track/#] <> " %"] &[N[maxDegree/100, 3]]
polCoefs = Join[{{{}, ∞}, {{}, 1}, {{}, 2}}, FoldList[(track = #2;
                polNext[#, #2]) &, {{{1, 9}}, 144}, Range[3, memorizeDegree]]];
(polMove[cur_, n_, m_] := If[# >= m, polCoefs[[m + 2]], Fold[polNext,
                             cur, Range[n + 1, m]]]) &[memorizeDegree]
polEval[{L_, den_}, odd_, x_] := If[den === ∞, 0, Fold[(First[#2] +
                                    x^2 Last[#2] #) &, 1, L] If[odd, x, 1]/den]

With[{stability = Reap[Fold[With[{polNext = polMove[#, #2 - 1, #2]}, track = #2;
      Sow[polEval[polNext, OddQ[#2], 1`10]]; polNext] &, {}, Range[0, maxDegree]];][[2, 1]]},
 polPrecs = Join[{∞, ∞}, Rest[extraPrecision - (Accuracy /@ stability - 10)]];
 maxAbsErr = SetPrecision[If[maxDegree >= 3, Times @@ Power[stability[[-4 ;;]],
                          {-1, 4, -6, 4}], {1/2, 5/72, 1/288}[[maxDegree + 1]]], 6];
 logSizes = -Log10[stability];];

The actual function:

fabius::limit = "Maximal degree reached; Probably failed to return " <>
                "the requested precision; The absolute error is at most `1`";
fabius[x_?(#1 ∈ Reals &), p_?(#1 ∈ Reals || # === ∞ &)] := Module[{refs, prec,
        cac1 = {{{1, 9}}, 144}, cac2 = 3, mod = Mod[SetPrecision[x, ∞], 2],
        iMax, pMax = Identity, xs, noDyadic, d2 = maxDegree + 2, sum = 0, sign,
        doN = MatchQ[x, _Rational | _Integer] && ByteCount[x] < 3000},
noDyadic = If[doN && IntegerQ[#], # > maxDegree, True] &[Log2[Denominator[mod]]];
refs = Position[IntegerDigits[BitXor[BitShiftRight[#, 1], #] &[FromDigits[
                  First[RealDigits[mod, 2, d2, 0]], 2]], 2, d2], 1][[All, 1]];
prec = p + logSizes[[Max[Min[Ceiling[1 - Log2[Min[mod, 2 - mod]]], Length[logSizes]], 1]]];
sum = SetAccuracy[sum, prec + extraPrecision]; prec += polPrecs[[refs]];
If[refs === {} && noDyadic, Message[fabius::limit, maxAbsErr]];
If[refs === {}, sum = 0, iMax = First[FirstPosition[prec, _?Negative, {All + 1}, {1}]] - 1;
 refs = refs[[;; iMax]]; prec = (prec[[;; iMax]] /. If[p =!= ∞, ∞ -> MachinePrecision, {}]);
 If[Not[doN], pMax = Function[u, SetPrecision[#, u] &][Max[prec]]];
 xs = FoldList[pMax[(1 - #) #2 - 1] &, pMax[Last[#] mod - 1], Most[#]] &[
                BitShiftLeft[1, Append[Differences[refs], First[refs] - 1]]];
 If[Not[doN], xs = MapThread[SetPrecision, {xs, prec}]];
 sign = {1, -1}[[FoldList[BitXor, ThueMorse[Floor[x/2]], BitAnd[Most[refs] - 1, 1]] + 1]];
 MapThread[(sum += #2 polEval[cac1 = polMove[cac1, cac2, cac2 = #], OddQ[#], #3]; 0) &,
  {refs - 2, sign, xs}]; If[noDyadic && iMax === All, Message[fabius::limit, maxAbsErr]];];
 (Remove[refs, prec, cac1, cac2, xs, pMax, mod, sign, sum, noDyadic, iMax, d2, doN];
  If[#2 && # != 0, SetPrecision[#, Min[Log10[Abs[#]] - Log10[maxAbsErr], Precision[#]]], #]
  ) &[sum, noDyadic]]

I'm not sure whether I should just remove the ancient code.

Answer 2

My version:

ClearAll[fabius];

fabius::usage = "fabius[x] gives the value of the Fabius function F(x) for a non-negative real argument x.";
fabius::nnreal = "Argument `1` should be a non-negative real number.";

Macros`SetArgumentCount[fabius, 1];
SyntaxInformation[fabius] = {"ArgumentsPattern" -> {_}};
Attributes[fabius] = {NumericFunction, Listable};

fabius[x : Underflow[] | Overflow[] | _DirectedInfinity] :=
(Message[fabius::nnreal, x]; Undefined);
fabius[x_] /; Im[x] != 0 || Re[x] < 0 := (Message[fabius::nnreal, x]; Undefined);

Module[{t, c, d, f},
  t[0] = 1; t[n_] := t[n] = (-1)^n t[Quotient[n, 2]];
  c[0] = 1; c[n_] := c[n] = Sum[c[n - m]/(2 m + 1)!, {m, n}]/(4^n - 1);
  d[n_] := d[n] = Sum[c[m]/(n - 2 m)!, {m, 0, n/2}]/2^(n (n + 1)/2);
  f[x_, p_] := f[x, p] = Module[{z = x, n, s, y = 0},
     While[z > 0,
       n = -Floor[Log2[z]]; z -= 2^-n;
       s = d[n] + Sum[d[n - m] 2^(m (m + 1)/2) z^m/m!, {m, n}]; y = s - y;
       If[p < Infinity && Abs[s/y] < 10^-p, Break[]]];
     SetPrecision[Abs[y], p]];

  fabius[x_] /; ExactNumberQ[x] && IntegerQ[Log2[Denominator[x]]] && 0 <= x <= 2 :=
  f[x, Infinity];
  fabius[x_] /; NumericQ[x] && x > 2 := fabius[Mod[x, 2]] t[Quotient[x, 2]];
  fabius[x_] /; NumericQ[x] && Im[x] == 0 && Re[x] >= 0 && Precision[x] < Infinity :=
  f[Rationalize[x, 0], Precision[x]];
  N[fabius[x_], {p_, _}] /; NumericQ[x] && Im[x] == 0 && Re[x] >= 0 := f[x, p];
  Derivative[n_][fabius][x_] := 2^(n (n + 1)/2) fabius[2^n x]];

Examples:

fabius[5/16]
(* 305857/2073600 *)

fabius[2/3]
(* fabius[2/3] *)

N[%, 60]
(* 0.819834885198518093044266564068975877132092199918258367435956 *)

fabius[0.775]
(* 0.953056 *)

Answer 3

Here is a slightly optimized version of Vladimir's code (which is an implementation of Rvachëv's method), based on results derived by Arias de Reyna in these two notes. (See also this related MO question.)

The two main ingredients in the following code are the optimized direct generation of the coefficients d in Vladimir's code (called ariasD[] in the implementation below), using Arias's recurrence, and the reformulation of the evaluation of the underlying approximating polynomials to Horner form. In addition, as I already noted in a comment above, SetPrecision[x, ∞] is more appropriate to use here than Rationalize[x, 0] for converting inexact arguments, since it directly generates the nearest dyadic rational (see this as well).

FabiusF::usage =
"FabiusF[x] gives the value of the Fabius function F(x) for a non-negative real argument x.";
Macros`SetArgumentCount[FabiusF, 1];
SyntaxInformation[FabiusF] = {"ArgumentsPattern" -> {_}};
SetAttributes[FabiusF, {NumericFunction, Listable}];

Derivative[n_Integer][FabiusF] := 2^(n (n + 1)/2) FabiusF[2^n #] &

(* https://mathematica.stackexchange.com/a/13245 *)
powerOfTwoQ[n_] := IntegerQ[n] && BitAnd[n, n - 1] == 0

FabiusF[Infinity] = Interval[{-1, 1}];
FabiusF[x_?NumberQ] /; If[0 <= Re[x] && Im[x] == 0,
                          powerOfTwoQ[Denominator[x]], 
                          Message[FabiusF::realnn, x]; False] := iFabiusF[x]

ariasD[0] = 1;
ariasD[n_Integer?Positive] := ariasD[n] = 
     Sum[2^((k (k - 1) - n (n - 1))/2) ariasD[k]/(n - k + 1)!, {k, 0, n - 1}]/(2^n - 1);

tri[x_] := Piecewise[{{2 - x, x > 1}}, x]

iFabiusF[x_] := Module[{prec = Precision[x], s = 1, y = 0, z = SetPrecision[x, Infinity],
                        n, p, q, tol, w},
       z = If[0 <= z <= 2, tri[z],
              q = Quotient[z, 2];
              (* can replace ThueMorse[] with the implementation
                 in https://mathematica.stackexchange.com/a/89351 *)
              If[ThueMorse[q] == 1, s = -1]; tri[z - 2 q]];
       tol = 10^(-prec);
       While[z > 0, n = -Floor[RealExponent[z, 2]]; p = 2^n; z -= 1/p; w = 1;
             Do[w = ariasD[m] + p z w/(n - m + 1); p /= 2, {m, n}];
             y = w - y;
             If[Abs[w] < Abs[y] tol, Break[]]];
       SetPrecision[s Abs[y], prec]]

Some examples:

FabiusF[Range[0, 10]]
   {0, 1, 0, -1, 0, -1, 0, 1, 0, -1, 0}

Table[{x, FabiusF[x]}, {x, {1/2, 1/4, 3/4, 1/8, 3/8, 5/8}}]
   {{1/2, 1/2}, {1/4, 5/72}, {3/4, 67/72}, {1/8, 1/288}, {3/8, 73/288}, {5/8, 215/288}}

{FabiusF[2/3], N[FabiusF[2/3], 30]}
   {FabiusF[2/3], 0.819834885198518093044266564069}

{FabiusF[π], N[FabiusF[π], 30]}
   {FabiusF[π], -0.993623885966442415174972734961}

and finally, a plot, just like what the OP asked for:

Plot[FabiusF[x], {x, 0, 24}, AspectRatio -> Automatic, Frame -> True]

As a bonus:

ParametricPlot[{FabiusF[t], FabiusF'[t]}, {t, 0, 4}]