(12C) Square Root

[Albert Chan]

Inspired by ∫(1/(a - cos(x)), x=0..pi) = pi/√(a²-1), I find another way to get 1/√(a²-1)

For x = 0 to pi, cos(x) = 1 to -1, we can approximate integral without cos(x) term:

∫(1/(a - cos(x)), x=0..pi) ≈ pi/a                     → 1/√(a²-1) ≈ 1/a

∫(1/(a - cos(x)), x=0..pi)
= ∫(1/(a - cos(x)) + 1/(a + cos(x)), x=0 .. pi/2);         // fold the integral
= 2*a * ∫(1/ (a² - cos(x)²), x = 0..pi/2)
= 4*a * ∫(1/ (2*a² - (1 + cos(2*x))), x = 0..pi/2);      // let y = 2*x
= 2*a * ∫(1/ (a₂-cos(y)), y = 0..pi);                           // let a₂=2*a²-1

Again, dropping cos(y), we have:

∫(1/(a - cos(x)), x=0..pi) ≈ 2*a/a₂ * pi          → 1/√(a²-1) ≈ (1/a)·(1+1/a₂)

Rinse and repeat, let a_k = 2*a_k-1² - 1, we have 1/√(a²-1) = (1/a)·(1+1/a₂)·(1+1/a₃) ...

Convergence is quadratic, example, with a=2

1/√3 = (1/2)·(1+1/7)·(1+1/97)·(1+1/18817)·(1+1/708158977) ...

Or, we flip both side:

√3 = 2·(1-1/8)·(1-1/98)·(1-1/18818)·(1-1/708158978) ...

function seq(a) -- sequence converging to sqrt(a*a-1), a > 1
local t = a
return function() a=2*a*a; t=t-t/a; a=a-1; return t end
end

lua> g = seq(2)
lua> for i=1,4 do print(i, g()) end
1 1.75
2 1.7321428571428572
3 1.7320508100147276
4 1.7320508075688774

This is equivalent to Newton's method for √N, N=a*a-1, guess x=a:

lua> N, x = 3, 2
lua> for i=1,4 do x=(x+N/x)/2; print(i, x) end
1 1.75
2 1.7321428571428572
3 1.7320508100147274
4 1.7320508075688772