Here's a solution using the sagetex package, which gives you access to a computer algebra system and Python programming.
\documentclass{article}
\usepackage{sagetex}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{pgfplots}
\pgfplotsset{compat=1.15}
\begin{document}
\begin{sagesilent}
def LSF(binexp):
a = .6666
L = [0]
U = [1]
M = [1]
for i in range(1,9):
M += [(1-a)*L[len(L)-1]+a*U[len(U)-1]]
if str(binexp)[i] == "1":
L += [M[len(M)-1]]
U += [U[len(U)-1]]
else:
U += [M[len(M)-1]]
L += [L[len(L)-1]]
return U[len(U)-1]
def BTD(mystr):
sum = 0
for i in range(1,9):
sum += int(mystr[i])*(.5)^i
return sum
xcoordsb = ['.{0:08b}'.format(i) for i in range(0,2^9)]
xcoords = [BTD(xcoordsb[i]) for i in range(0,2^9)]
ycoords = [LSF(xcoordsb[j]) for j in range(0,2^9)]
plotpoints = sorted([[xcoords[i],ycoords[i]] for i in range(0,2^9)], key=lambda k: [k[1], k[0]])
output = r""
output += r"\begin{tikzpicture}[scale=.7]"
output += r"\begin{axis}[xmin=0,xmax=1,ymin= 0,ymax=1,"
output += r"title={Lebesgue singular function, $a=.6666$}]"
output += r"\addplot[thin, blue, unbounded coords=jump] coordinates {"
for i in range(0,len(plotpoints)-1):
output += r"(%s,%s) "%(plotpoints[i][0],plotpoints[i][1])
output += r"};"
output += r"\end{axis}"
output += r"\end{tikzpicture}"
\end{sagesilent}
\sagestr{output}
\end{document}
The solution is shown running in Cocalc for when a=.6666. Changing the value
of a and then compiling will let you get other values.

My programming skills aren't strong, there's almost certainly an easier and more elegant ways to do this. I'm relying on an algorithm from Daniel Bernstein's thesis which I've linked to. The diagram on page 18 matches the output from Sage. The algorithm used is on page 19.
Sage is not part of the LaTeX distribution. The quickest way to get started using it is with a free Cocalc account. You can install Sage on your computer so you aren't relying on an internet account. Search this site for sagetex to get more information.