Yes, an acceleration is possible, a basic acceleration for all admissible scaling sequences, and even faster convergence the smoother the resulting wavelet scaling function is. This is based on knowing a related scaling function with closed formula for the function itself and also its Fourier transform, more specifically the B-splines. The recursion equation for the difference has a smaller contraction constant, down to about $2^{-A}$ if $A$ is the number of vanishing moments (db2 has $2$, db8 with sequence length $16$ has $A=8$).
My notation conventions
For my convenience, first switch notation to something more adapted to the task. Having a separate factor $\sqrt2$ is nice for exploring orthogonality of wavelets. For studying its convergence it is more helpful to have a factor $2$, which is also the reason to switch from $h$ to $m_0$ in the reference.
Shift operator Define $(\tau f)(x)=f(x-1)$. Then the scaling equation becomes $\phi(x)=2(a(\tau)\phi)(2x)$.
Dilation operator To remove the argument from the equation, additionally define a dilation or doubling operator $(Df)(x)=f(2x)$. Then the scaling equation becomes a fixed-point equation in function space $\phi=D(a(\tau)\phi)$. In a mash-up of these notations also use $\hat \tau(\xi)=\exp(-i\xi)$ for the basis shift.
Combining both operators results in $\tau D=D\tau^2$, and consequently the basic identity behind most of the following approach
$$
(1-\tau)D=D(1+\tau)(1-\tau).
$$
Fourier and Z transforms It makes sense to introduce the formal (Laurent) polynomial $a(z)$. Denote the Fourier factor of the convolution operator $a(\tau)$ as $\hat a(\xi)=a(\exp(-i\xi))$. In consequence $\hat\phi(\xi)=\hat a(\xi/2)\hat \phi(\xi/2)$.
Why slow convergence
The observed problem for the product is that the sequence $\xi/2^j$ converges fast toward $0$ and thus $\hat m_0(\xi/2^j)$ is ever closer to $1$, making for a very slow convergence.
Vanishing moments and factorization
We know that $a(1)=1$ from the integrability of the scaling function and $a(-1)=0$ as necessary condition for a continuous scaling function. This implies a factorization $a(z)=\frac{1+z}2p(z)$ with $p(1)=1$. Thus $p(z)=1+(1-z)q(z)$.
Comparing two scaling functions
Add a second pair $(\tilde a,\tilde\phi)$ of scaling sequence and function to the situation. Then with above conventions their difference can be represented as $\phi-\tilde\phi=(1-\tau)\eta$, and inserting the refinement equations
$$\begin{align}
(1-\tau)\eta &= 2D(a(\tau)\phi-\tilde a(\tau)\tilde\phi)\\
&=D(1+\tau)\Bigl(p(\tau)(\phi-\tilde\phi)+(p(\tau)-\tilde p(\tau))\tilde\phi)\Bigr)\\
&=(1-\tau)D\Bigl(p(\tau)\eta+(q(\tau)-\tilde q(\tau))\tilde\phi\Bigr)
\end{align}$$
This gives a recursion equation, one could say affine linear refinement equation, for $\eta$,
$$
\eta = D\Bigl(p(\tau)\eta+v(\tau)\tilde\phi\Bigr)
$$
with $v=q-\tilde q$.
In the Fourier picture this gives the fixed-point iteration
$$
\hat\eta(\xi)=\frac12\left(\hat p(\xi/2)\hat\eta(\xi/2)+\hat v(\xi/2)\hat{\tilde\phi}(\xi/2)\right)
$$
to be inserted in
$$
\hat\phi(\xi)=\hat{\tilde\phi}(\xi)+(1-\hat\tau(\xi))\hat \eta(\xi)\\
$$
Applying the Banach fixed-point theorem, this gives uniform convergence if $\sup|\hat p(\xi)|<2$.
This is now non-empty in the space of $p$ sequences, as the only a-priori restriction is $p(1)=1$.
Pointwise convergence is without restriction, as $p$ is polynomial with $\hat p(0)=p(1)=1$, so the sequence of factors $\hat p(\xi/2^n)/2$ converges to $1/2$.
A condition for convergence in $L^2(\Bbb R)$ is $\sup|\hat p(\xi)|<\sqrt 2$
B-splines
A class of known solutions to the refinement equation are the basis or B-splines $\beta_k$ with $\beta_0={\bf 1}_{[0,1]}$ the box function and recursively $\beta_{k+1}=\beta_0*\beta_k$. $\beta_1(x)=\max(0,1-|x-1|)$ is the hat function. The refinement sequence for $\beta_k$ is given via $a(z)=\left(\frac{1+z}{2}\right)^{k+1}$. The Fourier transform gets computed via
$$
\hat\beta_0(\xi)=\int_0^1e^{-ix\xi}\,dx=\frac{1-e^{-i\xi}}{i\xi}=e^{-i\xi/2}\frac{\sin(\xi/2)}{\xi}
\\
\implies\hat\beta_k=e^{-i(k+1)\xi/2}\sin^{k+1}(\xi/2)/\xi^{k+1}
$$
Examples with minimal support
Let $q$ be constant and take $\tilde\phi=\beta_1$, so that $\tilde p=\frac{1+z}2=1-\frac12(1-z)$. Then $v(z)=q_0+\frac12$. $|\hat p(\xi)|=\sqrt{(1+q_0)^2-2(1+q_0)q_0\cos(\xi)+q_0^2}$. The condition for uniform convergence is $|1+q_0|+|q_0|=\max(1,|1+2q_0|)<2$, that is, $q_0\in(-1.5,0.5)$.
The translation into an iteration for the Fourier transform follows directly.
The convergence is visually rapid in all cases, there is no visual difference between the 8th and (L=)50th iterate. How close $q_0$ is to $-0.5$ only determines how rapidly the amplitude falls, which is a measure for the smoothness of the scaling function itself.
The computation of a function table for a plot of the graph of the scaling function can be organized similarly.
As $p$ has degree 1, length 2, the iteration function $\eta$ has support $[0,1]$. The value $\eta(x)$ is used in the iteration formula for the values at $x/2$ and $(x+1)/2$.
$$
\pmatrix{\eta(x/2)\\\eta(x/1+1/2)}
=
\pmatrix{1+q_0\\-q_0}\eta(x)+v_0\pmatrix{x\\1-x}
\\
\pmatrix{\phi(x)\\\phi(x+1)}
=
\pmatrix{x\\1-x}+\pmatrix{1\\-1}\eta(x)
$$
This form shows that $\eta$ converges uniformly and thus to a continuous function for $\max(|1+q_0|,|q_0|)<1$, that is, for $q_0\in[-1,0]$. However, only for $q_0=-0.5$ is the limit Lipschitz continuous, moving away one can show Hölder continuity, with the index going to $0$ at the boundaries of the interval.
## start with the most basic subdivision
x,eta = np.array(([0,1],[0,0])
etas = [[x,eta]]
## iterate the eta recursion
for k in range(10):
eta = np.concatenate([(1+q)*eta+v*x, -q*eta+v*(1-x)])
x = np.concatenate([x/2,(1+x)/2])
etas.append([x,eta])
## construct and evaluate the scaling function
for k in [1,3,7]:
x,eta = etas[k]
phi = np.concatenate([x+eta, 1-x-eta])
x = np.concatenate([x,1+x])
plot(x,phi)
# etc
-O3 -ffast-math; now I get 300ns evaluation in float precision and 950ns evaluation in double. I imagine this can still be sped up . . . but it is no longer just totally insane. $\endgroup$