  We present algorithms to numerically evaluate Daubechies wavelets and scaling
functions to high relative accuracy. These algorithms refine the suggestion of
Daubechies and Lagarias to evaluate functions defined by two-scale difference
equations using splines; carefully choosing amongst a family of rapidly
convergent interpolators which effectively capture all the smoothness present
in the function and whose error term admits a small asymptotic constant. We are
also able to efficiently compute derivatives, though with a smoothness-induced
reduction in accuracy. An implementation is provided in the Boost Software
Library.
