# Compute the Fabius Function

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\$$ :

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

• 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.
– xnor
Sep 19, 2022 at 1:46
• This might be useful. Sep 19, 2022 at 2:04
• This formula might also help Sep 19, 2022 at 7:03

# Python, 94 bytes (@att)

lambda x:B.c[int(x*4**9)]/8**57
from numpy import*
B,*A=a=1,1.
exec('A*=2;B*=poly1d(A)*a;'*18)


Attempt This Online!

## Python, 99 bytes

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


Attempt This Online!

### Python, 133 bytes

lambda x:interp(x,r_[:2:len(B)+0j],B/max(B))
from numpy import*
A=a=poly1d([1,0])
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.

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

# PARI/GP, 56 bytes

-1 byte thanks to @att,

x->sum(i=0,x<<=20,(-1)^sumdigits(x\1-i,2)*i^20)/20!/4^95


Attempt This Online!

Returns a fraction.

\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.

• 2^190->4^95
– att
Sep 23, 2022 at 2:19

# Pyth, 41 39 bytes

J*hK20]ZV*^2KQ=JsM._XZJ^_1sjN2;ceJ^2sUK


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

{⊃(⌊⍵×8*7)↓(+\⊢,-)⍣22÷8*70}


Try it online!

As noted by CursorCoercer, the identity $$\f'(x)=2f(2x)\$$ holds when $$\0, 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]\$$.

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
approximation

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