10
\$\begingroup\$

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

Task

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
Improve this question
\$\endgroup\$
3
  • 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\$ Commented Sep 19, 2022 at 1:46
  • 3
    \$\begingroup\$ This might be useful. \$\endgroup\$ Commented Sep 19, 2022 at 2:04
  • 2
    \$\begingroup\$ This formula might also help \$\endgroup\$ Commented Sep 19, 2022 at 7:03

4 Answers 4

Reset to default
8
\$\begingroup\$

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!

-13 thanks @att

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.

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

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.

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.

Improve this answer
\$\endgroup\$
1
  • 1
    \$\begingroup\$ 2^190->4^95 \$\endgroup\$ Commented Sep 23, 2022 at 2:19
2
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
1
\$\begingroup\$

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<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
                           approximation

⊃(⌊⍵×8*7)↓          index into the samples
Improve this answer
\$\endgroup\$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.