How to calculate the double sine function via Sage to high precision?

edit

From a paper of Ruijsenaars, we can extract the following integral representation for the logarithm of the double sine function:

$$\log(S_2(w,z)) = \int_0^\infty \left( \frac{e^{-y(1+w-z)}-e^{-yz}}{(1-e^{-y})(1-e^{-wy})} - \frac{1}{y}(\frac{2z}{w}-\frac {1}{w} -1 ) \right) \frac{dy}{y} $$

In Sage, it is straightforward to implement this function using the procedure "numerical_integral". It works just fine and produces results in perfect alignment with the known values of the double sine function. (E.g., we have $S_2(w,1) = \sqrt{w}$).

My problem is, that I don't know how to increase the number of digits I can get using "numerical_integral". I only get around 16 digits in the decimal expansion, but for my applications I need at least 32.

Question 1: How to increase the number of digits in Sage for the output of the procedure "numerical_integral"? In other words, how can I increase the precision of the procedure "numerical_integral"?

Here is my code:

numerical_integral(( (e^(-(1+w-z)*y)-e^(-z*y))/((1-e^(-y))*(1-e^(-w*y))) -1/y*(2*z/w-1/w-1) )/y,0,+Infinity,max_points=200)

Thank you very much in advance for your help!

PS. Other suggestions to obtain high precision results are welcome as well of course...

PS2. I would be fantastic to have a proper implementation of the double sine function in Sage!