Estimating a quarter turn in radians

Compute two sequences $\begin{matrix}a_0&\dots&a_{n-1}\\b_0&\dots&b_{n-1}\end{matrix}$ of unsigned integers with
$$\begin{align}\begin{matrix}a_0&a_{n-1}\\b_0&b_{n-1}\end{matrix} &= \begin{matrix}1&0\\0&1,\end{matrix}\\\left|\begin{matrix}a_k&a_{k+1}\\b_k&b_{k+1}\end{matrix} \right| &= 1\;\text{for all $k.$}\end{align}$$Then $$\sum\frac1{a_ka_{k+1}+b_kb_{k+1}} \approx \frac\tau4.$$

The idea is as follows: Integrate 1 from 0 to $\frac\tau4.$ A small term $\Delta\theta$ can be approximated as $\sin \Delta\theta$. Take $\sin^2 \Delta\theta$ to be the spread between two adjacent vectors $\begin{matrix}a_k\\b_k\end{matrix}$, $\begin{matrix}a_{k+1}\\b_{k+1}\end{matrix}.$ This spread is $$\frac{(a_kb_{k+1}-a_{k+1}b_k)^2}{(a_k^2+b_k^2)(a_{k+1}^2+b_{k+1}^2)} = \frac1{(a_ka_{k+1}+b_kb_{k+1})^2+1} \approx \frac1{(a_ka_{k+1}+b_kb_{k+1})^2}.$$Thus $\Delta\theta\approx\frac1{a_ka_{k+1}+b_kb_{k+1}}.$

Here's an algorithm:

#include <iostream>
#include <iomanip>
using namespace std;

long double integrate(unsigned long a0, unsigned long b0,
                      unsigned long a1, unsigned long b1) {   
    if (a0*a1+b0*b1 > 100000)
        return 1.0L/(a0*a1+b0*b1);
    else
        return integrate(a0+a1, b0+b1, a1, b1) + integrate(a0, b0, a0+a1, b0+b1);
}

int main() {
    cout << "approx " << setprecision(20) << integrate(1, 0, 0, 1) << endl;
}

The output is 1.5707963268188070548, which equals $\frac\tau4$ to 9 decimal places.