$$\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.