Fabius.nb

(* Mathematica code for exact evaluation of the Fabius function and its derivatives at dyadic rationals *)
ClearAll[c, Fabius];

(* Coefficients in the Taylor—Maclaurin series of the Fourier image of the up(x) function *)
c[0] = 1;
c[n_Integer /; n > 0 && OddQ[n]] = 0;
c[n_Integer /; n > 0 && EvenQ[n]] := c[n] = 
   With[{k = n/2}, Sum[((-1)^(k - r) c[2 r])/(1 + 2 k - 2 r)!, {r, 0, k - 1}]/(4^k - 1)];

(* The Fabius function *)
Fabius[0] = 0;
Fabius[1] = 1;
Fabius[x_Rational /; 0 < x < 1 && IntegerQ[Log2[Denominator[x]]]] := Fabius[x] = 
   Module[{k = Numerator[x], n = Log2[Denominator[x]]}, 
    Sum[((-1)^(q + ThueMorse[p - 1])
       (1/2 - p + k)^(n - 2 q) c[2 q])/(4^q 2^((n - 1) n / 2) (n - 2 q)!), {p, 1, k}, {q, 0, n/2}]];
Fabius[x_ /; 1 < x <= 2] := Fabius[2 - x];
Fabius[x_ /; x > 2] := Fabius[Mod[x, 2]] (-1)^ThueMorse[Floor[x/2]];

(* The n-th derivative of the Fabius function *)
Derivative[n_Integer /; n > 0][Fabius][x_] := 2^(n (n + 1) / 2) Fabius[2^n x];