While trying to answer this question I used the series approach and obtained a recursive algorithm. While checking it numerically, I found it suffering from "catastrophic cancellation", i.e. loss of significant digits, and extremely slow convergence.
I would like to know if the algorithm can be improved or if I should just give it up as a bad idea (I was going to apply the same approach to other integrals).
With some simplifications from the previous answer, the algorithm is as follows:
$$0<k<1$$
$$ \begin{cases} T_0 = \dfrac{k^2}{\sqrt{1-k^2}} \\ P_0=1 \\ S_0=0 \end{cases}$$
$$ \begin{cases} T_n = \dfrac{1} {k^2} T_{n-1} -P_{n-1}\\ P_n=\left(1-\dfrac{1}{2n} \right)P_{n-1} \\ S_n=S_{n-1}+\dfrac{T_n}{2n-1} \end{cases} \quad,\quad n \geq 1$$
From another answer to the linked question, I know that:
$$\lim_{n \to \infty} S_n=\frac{k }{2 \sqrt{1-k^2}} \log \left(\frac{1+k}{1-k}\right)$$
But it had been extremely hard to get a good numerical agreement with the closed form, for the reasons stated above.
I understand why the cancellation happens: the expression for $T_n$ at each step contains subtraction of two numbers (apparently) close to each other in value.
This is what I get for $k=1/3$ in Mathematica, using $10^5$ steps of the algorithm and keeping $10^5$ digits for numerical evaluation at each step. Setting $10^4$ digits instead for the same number of steps makes the numerical result $0$.
In[804]:= Nm=100000;
k=1/3;
S0=0;
P0=1;
T0=k^2/Sqrt[1-k^2];
Do[T1=N[T0/k^2-P0,10^5];
P1=N[(1-1/(2 n)) P0,10^5];
S1=N[S0+T1/(2 n-1),10^2];
T0=T1;
P0=P1;
S0=S1,{n,1,Nm}];
N[S0,20]
N[(k Log[(1+k)/(1-k)])/(2 Sqrt[1-k^2]),20]
Out[810]= 0.12230925274427951206
Out[811]= 0.12253226793356839896
Note: I do not want closed forms for $T_n$ or $P_n$, I know them. I would like to know if the algorithm itself can be improved, without getting to the closed forms and then setting up some different algorithm based on them.