Programming Challenge: a classic trigonometry problem

[Albert Chan]

Perhaps better way, we solve for ratio p = v/u, then recover x.
This avoided worrying about hitting square root of negative numbers.

Again, assume a ≥ b > 0      → u ≥ v      → 0 < p ≤ 1

1/u + 1/v = 1/c      
u/c = (1 + 1/p)      .... (1)

(u²-v²) = (a²-b²)
u² * (1-p²) = (a²-b²)      .... (2)

(1) and (2): (1+1/p)² * (1-p²) = (a²-b²)/c² = k      .... (3)

Equation (3), p and c have positive relationship.

(c → min(c) = 0) ⇒ (p → 0)
(c → max(c) = 1/(1/a+1/b) = (a*b)/(a+b)) ⇒ (p → b/a)      .... (4)

We can estimate max(p), assuming (3) term (1-p²) = 1

From (3) and (4): p ≤ min(1/(√k - 1), b/a)

After solved for p, we recover x

From (2): u² = (a²-x²) = (a²-b²) / (1-p²)

x² = (b+p*a) * (b-p*a) / (1-p²)

lua> c = 15
lua> k = (a*a-b*b)/(c*c)
lua> P = fn'p: (1+1/p)^2*(1-p*p) - k'
lua> hi = min(1/(sqrt(k)-1), b/a)
lua> S.secant(P, 0.9*hi, hi, 1e-9, true)
0.75                               0.6936316463403318
0.6936316463403318      0.6922323943366336
0.6922323943366336      0.6923293702666201
0.6923293702666201      0.6923292025361991
0.6923292025361991      0.6923292025145776
0.6923292025145777
lua> p = _
lua> y = (b+p*a)*(b-p*a) / (1-p*p)
lua> sqrt(y) -- = x
15.987649008568583

Above example, b/a is used for hi estimate. For small c, it switched to the other.

lua> c = 1
lua> k = (a*a-b*b)/(c*c)
lua> hi = min(1/(sqrt(k)-1), b/a)
lua> S.secant(P, 0.9*hi, hi, 1e-9, true)
0.039281134636117175    0.03925422803226194
0.03925422803226194      0.03924967251136882
0.03924967251136882      0.039249677913927694
0.03924967791300026
lua> p = _
lua> y = (b+p*a)*(b-p*a) / (1-p*p)
lua> sqrt(y) -- = x
29.981993931474232