Ellipsoidal coordinates are an orthogonal coordinate system under which the Laplace equation can be solved by separation of variables. While this has many benefits over spherical coordinates for a variety of potential problems, computation in ellipsoidal coordinates is difficult. Most notably, high-order harmonics can lack closed-form solutions and the associated normalization constants require approximating a singular integral. We provide a method for computing normalization constants to machine precision using tanh-sinh quadrature which exhibits exponential convergence for a large class of functions with singular endpoints. Combined with previous efforts to make ellipsoidal harmonics more accessible, the result is a library which makes computation in ellipsoidal coordinates as accessible as computation in spherical coordinates. Finally, we apply our implementation to the mixed-dielectric solvation problem and provide work-precision analysis for the results.