Note that $f(x^2)=\frac12f(x)$, $x\in(0,1)$. Thus with $g(x)=1/f(e^{-x})$ we have $g(2x)=2g(x)$, $x\in(0,\infty)$. And with $h(x)=\ln g(e^x)-x$, $x\in\Bbb R$, we have $h(x+\ln 2)=h(x)$, which must therefore be periodic.
Numerically (summing $f$ only from $k=-50$ to $50$, but that should be pretty good when computing $h(x)$ for $0\le x\le\ln 2$), $h$ seems to make only a (sine-like?) oscillation between $\approx-0.3665228$ and $\approx-0.366503$ (both $\approx \ln\ln 2$). So if $h(x)=\ln\ln 2+\epsilon(x)$ with $\epsilon(x)\approx 0$, then $g(x)=e^{\ln\ln 2+\epsilon(\ln x)+\ln x}=v(x) x \ln 2$ with $\upsilon(x)=e^{\epsilon(\ln x)}\approx 1$, and then
$$f(x)=\frac1{g(-\ln x)}=\frac{-1}{v(-\ln x)\ln 2\ln x}.$$
The numerical results seem to indicate that $|\epsilon(x)|<9.885\cdot 10^{-6}$ so that $|\upsilon(x)-1|<9.885\cdot 10^{-6}$, i.e.,
$$f(x)\approx -\frac1{\ln 2\ln x} $$
with a relative error $<10^{-5}$. On the other hand, the oscillation does not seem to disappear when computing more summands of $f$, which means that this bound for the relative error
is pretty sharp.