The Fabius function is an example of a function that is infinitely differentiable everywhere, yet nowhere analytic.

One way to define the function is in terms of an infinite number of random variables. Specifically, given a sequence of independent random variables \$\{U_n\}\$, where each \$U_n\$ is uniform on the interval \$[0,2^{-n}]\$, the Fabius function \$f\$ is defined on \$[0,1]\$ as the cumulative distribution function of their sum \$\displaystyle\sum_{i=1}^\infty U_i\$ :

Plot of successive sums of the series of cdfs


Given a number \$x\in[0,1]\$, compute \$f(x)\$. Your answer should be accurate to within \$2^{-16}\approx0.00001526\$ on the given test cases.

Test cases

x       f(x)
0       0
0.163   0.012220888930835
0.3     0.129428260311965
0.5     0.5
0.618   0.733405658665583
0.75    0.930555555555556
1       1
  • 4
    \$\begingroup\$ Here's a restatement of the definition that I found helpful. Repeatedly apply the following operation starting with \$x\$: double it then subtract a uniformly random number from 0 to 1. Then, \$f(x)\$ is the chance it never goes negative. We can approximate this by applying the operation \$k\$ times and checking if the end result is non-negative, for some large \$k\$. Most of the time it quickly veers off towards infinity or negative infinity in a few steps. \$\endgroup\$
    – xnor
    Sep 19, 2022 at 1:46
  • 3
    \$\begingroup\$ This might be useful. \$\endgroup\$
    – alephalpha
    Sep 19, 2022 at 2:04
  • 2
    \$\begingroup\$ This formula might also help \$\endgroup\$
    – mousetail
    Sep 19, 2022 at 7:03

4 Answers 4


Python, 94 bytes (@att)

lambda x:B.c[int(x*4**9)]/8**57
from numpy import*

Attempt This Online!

Python, 99 bytes

lambda x:B.c[int(x*4**9)]/max(B.c) 
from numpy import*
for i in a*9:A*=2;B*=poly1d(A)*a

Attempt This Online!

-13 thanks @att

Python, 133 bytes

lambda x:interp(x,r_[:2:len(B)+0j],B/max(B))
from numpy import*
B=prod([((A:=A*A)/(a-1))[0]*(a+1)for i in[0]*17]).c

Attempt This Online!

Brute-force computation of the convolution of the given uniform densities. Uses polynomials to get convolution through the * symbol. Two little tricks: 1) x^n / (x-1) generates (discarding the residual 1) the "flat" polynomials we need 2) convolution with a constant is essentially prefix sum up to the middle, saving us a cumsum.

  • 1
    \$\begingroup\$ lambda x:B[int(x*2**17)]/max(B) \$\endgroup\$
    – att
    Sep 20, 2022 at 2:12
  • \$\begingroup\$ -4 bytes \$\endgroup\$
    – att
    Sep 22, 2022 at 3:31

PARI/GP, 56 bytes

-1 byte thanks to @att,


Attempt This Online!

Returns a fraction.

Based on this answer on Mathematics SE found by mousetail:

$$\begin{align} F(x)&=\lim_{n\to\infty} F_n(x)\\ F_n(x)&=\frac{2^{n(n+1)/2}}{n!}\,\sum_{y\in D_{n,x}}(-1)^{s(y)}\big(x-y\big)^n \end{align}$$

where \$D_{n,x}\$ is the set of all dyadic numbers in \$[0,x]\$ of the form \$m2^{-n}\$ for integers \$m\$, and \$s(y)\$ is the sum of the binary digits of \$y\$.

It seems that \$n=20\$ is sufficient for this challenge.

  • 1
    \$\begingroup\$ 2^190->4^95 \$\endgroup\$
    – att
    Sep 23, 2022 at 2:19

Pyth, 41 39 bytes


Try it online!

Note that in the try it online link, I have decreased the accuracy from 20 to 15 so that it finishes running in time. You'll have to take my word that I've tried it on my machine and the above code meets the given accuracy requirement.

This program essentially does N summations along the Thue–Morse sequence. This can be thought of as a sort of "discrete integral", in which case the fact that \$f'(x)=2f(2x)\$ tells us that each time we do this, we get our function back but with half the interval, half the magnitude, but twice the accuracy. We can easily deal with the magnitude at the end, and by starting with \$2^N\$ terms of the sequence, we end with the range from 0 to 1 with whatever accuracy we desire, No floating point or overflow errors occur along the way since all of this is done on python ints. \$N=20\$ is about right to hit the accuracy required.


APL (Dyalog Classic), 27 bytes


Try it online!

As noted by CursorCoercer, the identity \$f'(x)=2f(2x)\$ holds when \$0<x<\frac12\$, and can be used to characterize a unique extension of \$f\$ to \$[0,4]\$. Specifically, noting the symmetry of \$f\$ about \$x=0.5\$, we have \$f(x)=f(2-x)\$ for \$x\in[1,2]\$, then \$f(x)=-f(x-2)\$ when \$x\in[2,4]\$.

Plot of the extension of f from 0 to 4

From here, we can iterate \$\displaystyle f(x)=\int_0^{2x}f(t)\,dt\$ until we attain the desired precision:

          ÷8*70     start: 2^-210
                      one sample of an approximate-f on [0,2],
                      times a multiplicative factor.
(     )⍣22          repeat (22x) until precision goal reached:
   ⊢                  from n samples on [0,2],
    ,-                append its negation
                        => 2n samples on [0,4]
 +\                   cumulative sum
                        take the discrete integral.
                        the introduced multiplicative term of
                        1/∆x is cancelled by the start value.
                        => 2n samples on [0,2] of a better

⊃(⌊⍵×8*7)↓          index into the samples

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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