Programming Challenge: a classic trigonometry problem

[Gerson W. Barbosa]

(01-07-2014 06:03 AM)Thomas Klemm Wrote:  Wikipedia provides yet another solution: Crossed ladders problem

Yet another one:

'(15+d)/√((15+d)^2-225)=30/(√((20-4/3*d)^2-225)+√((15+d)^2-225))'

d: 2.72704740626

'E=5-4/3*d'

E: 1.36393679166

'x=√((15+E)^2-225)+√((15+d)^2-225)'

x: 15.9876490088

http://farm3.staticflickr.com/2860/11828...6eef_o.jpg

Not any more simple, but I just wanted to do it another way. Still using the solver, though.

Cheers,

Gerson.

[Werner]

If you take ladders of 3152m and 848m, and have them cross at a height of 585m, the distance between the walls is exactly 448m.
Or, take ladders of 4040m and 1160m, and have them cross at a height of 693m, then the distance between the walls is exactly 800m.
;-))
Werner

[Albert Chan]

(01-01-2014 08:44 AM)Thomas Klemm Wrote:  I used additional variables:
\[
u=\sqrt{a^2-x^2} \\
v=\sqrt{b^2-x^2}
\]...
\[
\frac{1}{u}+\frac{1}{v}=\frac{1}{c}
\]

Avoiding infinities, when u or v is zero, let \(f = c\,(u+v) - uv\)

\( \displaystyle \frac{d}{dx}\sqrt{k-x^2} = \frac{-x}{\sqrt{k-x^2}}\)

For both u and v, if we increase x, slope numerator goes bigger, but denominator goes smaller.
To straighten f curve, let \(y = x^2\)

\( \displaystyle \frac{d}{dy}\sqrt{k-y} = \frac{-1/2}{\sqrt{k-y}}\)

Assume a ≥ b > 0, and 1/c > 1/a + 1/b

\(f(y = 0) = c\,(a+b) - ab = abc \left(\frac{1}{a} + \frac{1}{b} - \frac{1}{c} \right) < 0\)

\(f(y = b^2) = c\,(u + 0) - u×0 > 0\)

We have a guaranteed real root between the y extremes.
Solved with Ridder's method (FYI, Free42/Plus42 SOLVE also use Ridder's method)

Note: each updated interval required cost of 2 function calls.
Note: last result is secant root of last interval, fall back to bisection, or f(x)=0

lua> S = require'solver' -- translated from Ridder solver Python code
lua> f = fn'y: local u,v = sqrt(a*a-y), sqrt(b*b-y); c*(u+v)-u*v'
lua> a, b, c = 40, 30, 15
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
0                      268.5639964554787
134.28199822773934     255.620329669562
194.95116394865067     255.6049251825406
225.27804456559562     255.60492082148892
255.6049208211841
lua> sqrt(_) -- = x
15.987649008568587

This work for full range c = 0 .. a*b/(a+b)

lua> c = 1
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
885.4877793349734      900
897.8097220135339      900
898.9165749684195      900
898.9165749684195      898.9214889155777
898.9199601065723      898.9214889155777
898.9199601065723      898.9199601069578
898.9199601069578
lua> sqrt(_) -- = x
29.981993931474232

lua> c = 17
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
0                      19.936187659463883
9.968093829731941      18.368152011527595
14.168122920629768     18.368139023904078
18.368139023879372
lua> sqrt(_) -- = x
4.285806694646805

---

Ridder's method, reduced interval may not translate to less iterations.
Ridder's method convergence rate may degenerate to bisection.

1/v = 1/c - 1/u ≤ 1/c - 1/a
v ≥ (a*c/(a-c))
(v² = b²-y) ≥ (a*c/(a-c))²

y ≤ (hi = b² - (a*c/(a-c))²)

max(c) = max(u*v/(u+v)) = a*b/(a+b)
min(c) = min(u*v/(u+v)) = u*0/(u+0) = 0

When y is small, (|u'|, |v'|) are relatively small, due to dominated constants, (a,b).
When y close to b², |v'| goes infinite, and c goes to 0 like fallen down a cliff.
This explains why f(y) look like a straight line, except when y very close to b²

We expect u*v/(u+v) concave-down, shaped like a quarter ellipse.

c/max(c) ≥ 1 - y/b²      // equality at the edges, y = 0 or b²

y ≥ (lo = b² * (1 - c*(a+b)/(a*b)))

Previously we simply use y = 0 .. b² = 900

lua> c = 15
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
112.5                  324
218.25                 255.6187346135275
236.93436730676376     255.6049211845503
246.26964424565705     255.60492082118668
250.93728253342186     255.6049208211843
255.6049208211839      255.6049208211843
255.60492082118404

lua> c = 1
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
847.5                  898.9480604865221
898.8848145203102      898.9480604865221
898.9199598288078      898.9480604865221
898.9199598288078      898.9199601069668
898.9199599678873      898.9199601069578
898.9199601069577      898.9199601069578
898.9199601069578

lua> c = 17
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
7.499999999999973      25.89792060491493
16.69896030245745      18.368146413229564
17.53355335784351      18.368139023879934
17.95084619086172      18.368139023879515
18.368139023879092     18.368139023879515
18.3681390238793

---

Reduced interval may not worth the trouble, unless we use it for secant's method.
I would ues 2 maximum guess, and iterate down, since f(y<0) is valid.

lua> S.secant(f, b*b, b*b-(a*c/(a-c))^2, 1e-9, true)

a/(a-c) = 1/(1-(c/a)) = 1 + (c/a) + (c/a)² + ... ≥ 1
We could use a simpler, but worse maximum guess.

lua> S.secant(f, b*b, b*b-c*c, 1e-9, true)

Or, we make f(y>b*b) vaild too, as long as we don't change the root.

lua> ssqrt = fn'x: x<0 and -sqrt(-x) or sqrt(x)'
lua> f = fn'y: local u,v = ssqrt(a*a-y), ssqrt(b*b-y); c*(u+v)-u*v'

[Massimo Gnerucci]

[rprosperi]

(03-07-2023 07:01 PM)Massimo Gnerucci Wrote:  

Well... 9 anyhow.

But Alvin Lee (of TYA) still has the record for fastest guitar solo in "I'm Going Home" and it's been nearly 54 years. Thanks for the refresh...

[cortopar]

I’d love to change the world…

[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

[Albert Chan]

(01-07-2014 06:03 AM)Thomas Klemm Wrote:  Wikipedia provides yet another solution: Crossed ladders problem
\[ x^3(x-c)=1 \]
Which leads to:
\[ x=\frac{1}{x^3}+c \]

Note: X ≈ C if C is big; X ≈ 1 if C is tiny.      → X = [C, C+1]

lua> a, b, c = 40, 30, 15
lua> D = sqrt(a*a-b*b)
lua> C = 4*c / D
lua> S = require'solver'
lua> S.secant(fn'X: (X-C)-X^-3', C, C+1, 1e-9, true)
3.2677868380553634      2.3488985640799114
2.3488985640799114      2.3451475080781123
2.3451475080781123      2.345304807136498
2.345304807136498        2.345304763754943
2.345304763754943        2.345304763754418
2.345304763754418
lua> X = _
lua> u = D*(X+1/X)/2      -- FYI: u+v = D*X; u-v = D/X
lua> sqrt((a+u)*(a-u))     -- = x
15.987649008568585

[Albert Chan]

(03-09-2023 04:19 AM)Albert Chan Wrote:  Note: X ≈ C if C is big; X ≈ 1 if C is tiny.      → X = [C, C+1]

X guesses have numerical issues if C = ε.

F(X) = (X-C) - 1/X^3

X1 = ε      → Y1 = 0 - 1/ε^3, which is huge
X2 = 1+ε  → Y2 = 1 - 1/(1+ε)^3 ≈ 0, numerically has order of machine epsilon.

Secant line slope = (Y2-Y1)/(X2-X1) ≈ (0-(-1/ε^3)) / ((1+ε)-ε) = 1/ε^3

newX = X1 - Y1/slope ≈ ε - (-1/ε^3) / (1/ε^3) = 1+ε = X2

It is unable to improve. Worse, it assumed convergence, wrongly report X2 as the answer.

lua> C = 1e-6
lua> S.secant(fn'X: (X-C)-X^-3', C, C+1, 1e-9, true) -- BAD!!!
1.000001      1.000001
1.000001

We wanted (X1, X2) closer, so that (Y1, Y2) about same magnitude.
Here is my solution, X1, X2 = C+1, C+exp(-C)

If C is tiny, C + exp(-C) ≈ 1 + C²/2

This make (X1,X2) gap small, exactly what we wanted.

lua> C = 1e-6
lua> S.secant(fn'X: (X-C)-X^-3', C+1, C+exp(-C), 1e-9, true) -- GOOD
1.0000000000005        1.000000250000375
1.000000250000375      1.0000002500000937
1.0000002500000937

If C is huge, C + exp(-C) ≈ C (slightly bigger than C, likely a better guess)

Redo previous post example, for X

lua> a, b, c = 40, 30, 15
lua> D = sqrt(a*a-b*b)
lua> C = 4*c / D
lua> S.secant(fn'X: (X-C)-X^-3', C+1, C+exp(-C), 1e-9, true)
2.37132791792656       2.3441836383192522
2.3441836383192522     2.345306973289791
2.345306973289791      2.3453047639451214
2.3453047639451214     2.345304763754418
2.345304763754418

[Gil]

Just tried today to tackle the problem.

Just wrote the equations of the red line and green line.

Had them cross at y=15

From there, I could eliminate the intercept of green equation.

And used the solver of the HP50G for the found Polynomial expression.

But I get two real values (solutions ?):
15.98 and 28.3565769436.

[Gil]

Yred = slope × "values on X-ax)
& slope red = opp/adj = sqrt (40² - x²) / x
Yred (at Xo, corresponding to yo = 15) <==>
sqrt (40² - x²) / x × Xo = 15
—>Xo= 15 × x/sqrt (40² - x²) (eq 1)

Green (at Xo):
Intercept + slope × Xo = 15 (eq 2)
& slope Red =- sqrt (30²-x²) / x
& intercept = sqrt (30²-x²)

EQ 2 Green (at Xo): <==>

sqrt (30²-x²) - [sqrt (30²-x²) / x] × Xo = 15
sqrt (30²-x²) - [sqrt (30²-x²) / x] × x/sqrt (40² - x²) = 15 (with eq 1)

Solve, getting rid of sqrt and writing x² = X

And we find:
'X^4.+-4100.*X^3.+5755000.*X^2.+-3091500000.*X+478406250000.'

Then:
[ 1. -4100. 5755000. -3091500000. 478406250000. ] in HP50G polynomial solver:

—>Roots:
[ (255.604920821,0.) (804.095455958,0.) (1520.14981161,129.64234902) (1520.14981161,-129.64234902) ]

Question :
Is sqrt of 804.095455958 not an acceptable solution?

[Albert Chan]

(03-09-2023 04:57 PM)Gil Wrote:  But I get two real values (solutions ?): 15.98 and 28.3565769436.

The bigger width corresponds to ratio p = v/u < 0
The ladders don't actually cross, but its line extensions does.

(03-08-2023 04:09 PM)Albert Chan Wrote:  (1) and (2): (1+1/p)² * (1-p²) = (a²-b²)/c² = k      .... (3)

We can estimate location of negative p

(1/(-p) - 1)² * (1-p²) = k
(1/(-p) - 1) ≥ √k

(-p) ≤ 1/(√k + 1)

lua> a, b, c = 40, 30, 15
lua> k = (a*a-b*b)/(c*c)
lua> P = fn'p: (1+1/p)^2*(1-p*p) - k'
lua> hi = -1/(sqrt(k)+1) -- guess for p (with negative sign)
lua> S = require'solver'
lua> S.secant(P, 0.9*hi, hi, 1e-9, true)
-0.3618162034940813        -0.34877707736298674
-0.34877707736298674      -0.34699739631288035
-0.34699739631288035      -0.3471286984646765
-0.3471286984646765        -0.34712756456695126
-0.34712756456695126      -0.3471275637878545
-0.34712756378785914
lua> p = _
lua> v = c * (1+p)
lua> u = v / p
lua> u, v
-28.21178023525376      9.793086543182113

Note that u < 0, with ladder *below* ground level.

---

Let xu = width from ladder "u" to where ladders cross
Let xv = width from ladder "v" to where ladders cross

x = xu + xv

From similar triangles, we have:

p = v / u = xv / xu

Above example, width numbers:

lua> x = sqrt((a+u)*(a-u)) * (u<0 and -1 or 1)
lua> xu = x / (1+p)
lua> xv = p * xu
lua> x, xu, xv
-28.356576943590095     -43.43356430868842      15.07698736509832

Extended lines cross 15 length units from ladder "v", in the opposite direction.

Update:

For extended-line crossed solution, I updated sign of x matching sign of u.
It is more apparent now that ladders never crossed.
It is also more symmetric, with sign of xu also matching sign of u.

Angles between ladders formula is simpler, without needing abs()
Let γ = angles between ladders, or its line extensions.

Area of quadrilateral = ½ a b sin(γ) = ½ (u+v) x = ½ (u v / c) x

sin(γ) = (u v x) / (a b c)

lua> (u*v*x) / (a*b*c) -- sin(γ)
0.43524258739970595
lua> deg(asin(_)) -- γ, in degrees
25.80073109384667

[Gil]

Thanks for your detailed answer, Albert.

My post illustrates the danger of a dummy "knowing" something, unduly risking "solutions", vs knowing nothing and standing wisely apart.

Regards,
Gil

[Pekis]

Interesting enough, on the crossed ladders Wikipedia page, in order to fold a sheet of paper in thirds if there is already a mark on the half height:

Starting from 1/c = 1/sqrt(a^2 - d^2) + 1/sqrt(b^2 - d^2)
(where d is the distance between the ladders)

Knowing (with Pythagore):
a^2 = d^2 + 1^2 (where 1 is the relative sheet height) => a^2 - d^2 = 1
b^2 = d^2 + (1/2)^2 (where 1/2 is the relative sheet half height) => b^2 -d^2 = 1/4

So
1/c = 1/sqrt(1) + 1/sqrtr(1/4) = 3 => c = 1/3
i.e the intersection point will be at the the third of the relative sheet height, for any relative sheet width !

[Albert Chan]

(03-08-2023 04:09 PM)Albert Chan Wrote:  (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)

With maximum guess, it is possible secant line over-shoot passed y-axis.
It may blow up, or converge to the "wrong" root.

It would be better if we have a minimum guess instead.
This way, iterations converge away from y-axis, making it more robust.

Let q = 1/p = u/v      // p max guess = q min guess

Q = (1+q)^2*(1-q^-2) - k

Q = 0       ⇒ |q| and k have positive relationship.

k = 16 ε³  ⇒ revert taylor series, q ≈ -horner([2,2,1],ε), +horner([2,2,1],ε³)
k → ∞      ⇒ Q ≈ q^2 - k = 0,     q ≈ ±√k

I have modified secant verbose option to display iterations, 1 per line.

I also added option for negative eps, to swap (x1,x2) if |f(x1)| < |f(x2)|
"Worse" point is moved as the first, thrown-away ASAP.

lua> S = require'solver'
lua> a, b, c = 40, 30, 15
lua> k = (a*a-b*b)/(c*c)
lua> Q = fn'q: (1+q)^2*(1-q^-2) - k'
lua> eps = -1e-9 -- swap worse point as first

lua> lo = max(sqrt(k)-1, a/b) -- ladder crossed solution
lua> S.secant(Q, lo, 1.1*lo, eps, true)
1.3333333333333333 -- = lo
1.4666666666666668
1.4445731516722746
1.444399371124978
1.444399566520114
1.4443995665182763

lua> q = _
lua> u = c*(1+q)
lua> x = sqrt((a-u)*(a+u))
lua> x, u, a -- right triangle
15.987649008568585      36.665993497774146      40

Extended-line crossed solution (negative q), we may just flip the sign of lo, and solve.
It take more iterations than with good guess, but it can be done. Over-shoot is harmless.

lua> lo = -lo -- bad guess, but OK
lua> S.secant(Q, lo, 1.1*lo, eps, true)
-1.3333333333333333 -- = lo
-1.4666666666666668
-7.344688904660508
-1.9135094120851166
-2.263168070942035
-3.2061293483685995
-2.8172286303461243
-2.875248024422489
-2.8808890380615035
-2.8807851830598548
-2.8807853490119513
-2.880785349016917
-2.880785349016917

lua> lo = -max(sqrt(k)+1, a/b) -- good guess
lua> S.secant(Q, lo, 1.1*lo, eps, true)
-3.040217628114033
-2.7638342073763935 -- = lo, swapped here
-2.875471300898932
-2.880971458358896
-2.880785063180847
-2.8807853490015667
-2.880785349016917
-2.880785349016917

lua> q = _
lua> u = c*(1+q)
lua> x = sqrt((a-u)*(a+u))
lua> x, u, a -- right triangle
28.3565769435901        -28.211780235253755      40

solve (1+q)^2 * (1-q^-2) * UnitStep[q^2-1] = 28/9, for q

Comment: since k ≥ 0, plot only show |q| = |u/v| ≥ a/b ≥ 1

---

From my next post, I suggested a correction, inside sqrt: (v/c)² = (1+1/q)²

Correction = 1 if q=±∞, but k ≈ q², even bigger. Correction doesn't matter.
Correction = 4 if q=1, it is worthwhile to use. Bonus: now hi ≥ 1

lua> hi = -1 + sqrt(k+4)
lua> S.secant(Q, hi, hi*0.9, eps, true)
1.6666666666666665 -- = hi
1.5
1.4442155853072434
1.4444000308015605
1.4443995665229092
1.4443995665182763
1.4443995665182763

Correction = 0 if q=-1., lo = -1 - sqrt(k) seems OK.
But when k is tiny, lo guess might be "too good".
Q'(q=-1) = 0. Iterations may bounce around wildly.
Last 2 points are extrapolated again for secant root, which may be garbage.
Going from hi side, steadily reach the true root is safer. (Q is concave-up)

lua> hi = -1 - sqrt(k+1)
lua> S.secant(Q, hi, hi*0.9, eps, true)
-3.0275875100994063 -- = hi
-2.7248287590894655
-2.874159866479985
-2.881098675537131
-2.880784748823651
-2.8807853489626534
-2.880785349016917
-2.880785349016917

[Albert Chan]

Here is an interesting symmetric pair of q guesses, derived at the same time.

1/v = 1/c - 1/u

Remember v is positive, but u can be both (u has sign of q): edge(u) = ±a
Let z = c/a

edge(1/v) = 1/c ± 1/a = (1 ± c/a) / c
edge(v/c) = (1 ± z)^-1

u = ±√(a² - y) = ±√(a² - (b² - v²)) = c * (1 + q)

(1 + q) = ±√((a²-b²)/c² + (v/c)²)

edge(q) = -(1 ± √(k + (1±z)^-2))

We redo previous post example, solve for both Q real roots.
I named negative q guess as "hi", only because it is bigger in size.

lua> S = require'solver'
lua> a, b, c = 40, 30, 15
lua> k = (a*a-b*b)/(c*c)
lua> Q = fn'q: (1+q)^2*(1-q^-2) - k'
lua> eps = -1e-9 -- swap worse point as first

lua> z = c/a
lua> lo = -(1 - sqrt(k + (1-z)^-2)) -- positive q guess
lua> hi = -(1 + sqrt(k + (1+z)^-2)) -- negative q guess

lua> S.secant(Q, lo, lo*1.01, eps, true)
1.3814094799322336 -- = lo
1.395223574731556
1.4441708057247238
1.4443988616014456
1.4443995665095293
1.4443995665182763

lua> S.secant(Q, hi*0.99, hi, eps, true)
-2.907888028932993 -- = hi
-2.878809148643663
-2.8807700072753484
-2.8807853577712663
-2.8807853490168784
-2.880785349016917
-2.880785349016917

[Albert Chan]

(03-10-2023 03:53 PM)Albert Chan Wrote:  Correction = 0 if q=-1., lo = -1 - sqrt(k) seems OK.
But when k is tiny, lo guess might be "too good".
Q'(q=-1) = 0. Iterations may bounce around wildly.

With Secant's method, we want to avoid zero slope situation.
And, if q ≈ -1, u = c*(1+q) suffer subtraction cancellation.

For Q negative root, perhaps we should let r = -(1+q) ≥ 0

(1+q)^2 * (1 - q^-2) = k
r^2 - r^2/(1+r)^2 = k

Let s = r/(1+r) ≤ min(r, 1)

(r-s) * (r+s) = k
(r-s) = r - r/(1+r) = r*s = k/(r+s)

r*r ≥ r*s = k/(r+s) ≥ k/(r+r)
2 r³ ≥ k
r ≥ ³√(k/2)

From previous post, we have another: r ≤ √(k+1)      (*)

Solve for r is more robust (no zero slope issue), faster (less iterations), and better accuracy.
This setup is stable. We can even go for full convergence!

lua> S = require'solver'
lua> eps = 0 -- full convergence
lua> R = fn'r,s: s=r/(1+r); r*s - k/(r+s)'

lua> k = 1e9
lua> r = S.secant(R, cbrt(k/2), sqrt(k+1), -eps, true)
793.7005259840997
31622.776617495183
31622.776617495136
31622.776617494186
31622.776617494186

lua> k = 1e-9
lua> r = S.secant(R, cbrt(k/2), sqrt(k+1), -eps, true)
0.0007937005259840997
1.0000000005
-- points swapped
0.0007937020238029177
0.0007940155478145599
0.0007940155478967004

For k = 1e-9 --> R = 5 calls

Solve Q(k = 1e-9) = 0 for negative q, for comparison:
hi = -1-sqrt(k+1), 0.9*hi .. hi --> Q = 32 calls
lo = -1-sqrt(k),     lo .. 1.1*lo  --> Q = 795 calls!

(*) we also have r ≥ √k, so technically r ≥ max(√k, ³√(k/2))
But, I don't want too complicated limit, especially since calculating R is so cheap.
We are only interested in asymptote estimate when k in tiny, thus cube root only.

(01-07-2014 06:03 AM)Thomas Klemm Wrote:  Wikipedia provides yet another solution: Crossed ladders problem
\[ x^3(x-c)=1 \]

Another way is to solve above, for negative X root.
Let Y = -1/X > 0, we have:

Y^4 = C*Y + 1      → Y = ⁴√(C*Y+1)

If C is big, we can drop +1, and get Y ≈ ³√C
The formula is so stable, that same guess is good enough for small C
We can even use it as iteration formula, without using solver!

Again, we solve for r = -(1+q), indirectly from y, also full convergence.
To simplify, we assume c=1. In other words, all variables are dimensionless.

lua> Y = fn'y: y - (C*y+1)^0.25'

lua> k = 1e9
lua> C = 4/sqrt(k); lo = cbrt(C)
lua> y = S.secant(Y, lo, lo*1.01, -eps, true)
0.05019802884366822
0.0507000091321049
1.000031623628705
1.0000316222766017
lua> (y+1/y) * (2/C) -- = r
31622.776617494186

lua> k = 1e-9
lua> C = 4/sqrt(k); lo = cbrt(C)
lua> y = S.secant(Y, lo, lo*1.01, -eps, true)
50.19802884366822
50.7000091321049
-- points swapped
50.198031475628454
50.19803147889932
50.198031478899324
lua> (y+1/y) * (2/C) -- = r
0.0007940155478967004

[Albert Chan]

(03-12-2023 05:18 AM)Albert Chan Wrote:  edge(1/v) = 1/c ± 1/a = (1 ± c/a) / c

There are many ways to look at effect at the edges.
edge(v/c) = (1 ± c/a)^-1, with -1 power, gives minimum guess for q

Here is another way.

edge(v/c) = edge(1+1/q) = (1 ± b/a), should gives maximum guess for q

(1+q) = ±√((a²-b²)/c² + (v/c)²), with both ways for edge(v/c), we have:

lua> S = require2'solver'
lua> Q = fn'q: (1+q)^2*(1-q^-2) - k'
lua> a, b, c = 40, 30, 15
lua> k = (a+b)*(a-b)/(c*c)
lua> eps = -1e-9 -- swap worse point as first

lua> S.secant(Q, -1+sqrt(k+(1-c/a)^-2), -1+sqrt(k+(1+b/a)^2), eps, true)
1.3814094799322336
1.4846752526459288
1.4445475888872654
1.4443992821167904
1.4443995665205573
1.4443995665182763
1.4443995665182763

lua> S.secant(Q, -1-sqrt(k+(1+c/a)^-2), -1-sqrt(k+(1-b/a)^2), eps, true)
-2.907888028932993
-2.781463193869329
-- points swapped
-2.8799914833475135
-2.8807791881861102
-2.8807853504286443
-2.8807853490169144
-2.880785349016917
-2.880785349016917

The 2 guesses are asymptote of extreme cases.
Minimum q guess is better when k is huge
Maximum q guess is better when k is tiny.

For middle ground, we may use geometric mean of 2 edge(v/c):

q ≈ -1 ± sqrt(k + (±a+b)/(±a-c))

[Albert Chan]

(03-07-2023 06:43 PM)Albert Chan Wrote:  max(c) = max(u*v/(u+v)) = a*b/(a+b)
min(c) = min(u*v/(u+v)) = u*0/(u+0) = 0

When y is small, (|u'|, |v'|) are relatively small, due to dominated constants, (a,b).
When y close to b², |v'| goes infinite, and c goes to 0 like fallen down a cliff.
This explains why f(y) look like a straight line, except when y very close to b²

We expect u*v/(u+v) concave-down, shaped like a quarter ellipse.

c/max(c) ≥ 1 - y/b²      // equality at the edges, y = 0 or b²

y ≥ (lo = b² * (1 - c*(a+b)/(a*b)))

Something magical just happened

From plots, it is obvious w = u*v/(u+v), as a function of y, concave down.
(location where w = c is the y we seek: 1/u + 1/v = 1/c)

Out of curiosity, I tried to plot w^2 vs y.
It look like a striaght line, only noticeably concave down when y close to b².
If true, we have an even better low y:

(c/max(c))² ≥ 1 - y/b²      → y ≥ (lo = b² * (1 - (c/max(c))²))

I tried to prove this, (w^2)'' < 0, but gave up because too many variables to juggle.
We can make length dimensionless, say all as ratio to a, but that still leave b and c.

My previous post, by plain luck, solved this! (all variables combined as ratio, q = u/v)
To make prove robust, we need to fill in inequality sign ... but thats it!

All is needed is to confirm both ways (low y from q, or concave down w²) are equivalent.

CAS> assume(a>=b); assume(b>=c); assume(c>=0)
CAS> k := (a+b)*(a-b)/(c*c)
CAS> cmax := a*b/(a+b)
CAS> u := c*(1+q)
CAS> yq := a*a - u*u

CAS> q := -1 + sqrt(k + (1+b/a)^2)  /* q high guess, from previous post */
CAS> y := b*b * (1 - (c/cmax)^2)      /* lo = b*b - c*c * (1+b/a)^2 */
CAS> simplify(y - yq)                  → 0

These are the y equivalent of previous post other q guesses.

CAS> q := -1 + sqrt(k + (1-c/a)^-2)  /* q low guess, from previous post */
CAS> y := b*b - (a*c/(a-c))^2           /* hi = b*b - c*c * (1-c/a)^-2 */
CAS> simplify(y - yq)                  → 0

CAS> q := -1 + sqrt(k + (a+b)/(a-c)) /* q middle guess, from previous post */
CAS> y := b*b - c*c * (a+b)/(a-c)
CAS> simplify(y - yq)                  → 0

Last y expression is simple, and pretty good rough guess.

CAS> y(a=40., b=30., c=15.)      → 270.0
CAS> sqrt(Ans)                          → 16.4316767252      // true x = 15.9876490086

[Albert Chan]

We could solve for s = (v/c), then v = c*s, then y = b^2 - v^2

c/u + c/v = c/c

±1/√(k+s²) + 1/s = 1       // ± to cover extend-line crossed solution
±s/√(k+s²) + 1 = s          // scaled to make it more linear     
|s-1| - s/√(k+s²) = 0       // formula V shaped, very straight for s≥1    (*)

This setup does not have secant's root overshoot problem.
It will not be hit with square root of negative number.

q ≤ -1      → s = (v/c) = (1+1/q) = 0 .. 1
q ≥ 1       → s = (v/c) = (1+1/q) = 1 .. 2

We can simply use guess [1±0.3, 1±0.7], instead of [(1-c/±a)^-1, (1+b/±a)]

Formula is concave up (2nd derivative > 0), it will not accidentally solve the "wrong" root.

CAS> factor(diff(abs(s-1) - s/sqrt(k+s*s), s, 2))      → \(\displaystyle \frac{3\,k\,s}{(k+s^2)^{5 \over 2}}\)      (**)

lua> a,b,c = 40,30,15
lua> k = (a+b)*(a-b)/(c*c)
lua> S1 = fn's: abs(s-1) - s/sqrt(k+s*s)'
lua> eps = -1e-9 -- swap worse point as first
lua> S = require'solver'
lua> S.secant(S1, 1.3, 1.7, eps, true)
1.3
1.7
1.6919260149527686
1.692328847882308
1.692329202531022
1.6923292025145777
1.6923292025145775
lua> v = c * _
lua> y = (b+v) * (b-v)
lua> y, sqrt(y) -- ladders crossed width
255.6049208211842      15.98764900856859

lua> S.secant(S1, 0.3, 0.7, eps, true)
0.3
0.7
0.6541746293109201
0.6528669490148099
0.6528724368426823
0.6528724362121412
0.6528724362121409
lua> v = c * _
lua> y = (b+v) * (b-v)
lua> y, sqrt(y) -- extended line crossed width
804.0954559577455      28.356576943590095

(*) s ≤ 1 side is not as straight, especially around s ≈ 0.
To reduce iterations, for s < 1 root with tiny k, use 3 point quadratic fit solver.

(**) there is a discontinuity at s=1. It should have another term, +abs(s-1)''
For our purpose of checking concave behavior, this subtlety can be ignored.

(03-17-2023 03:12 PM)Albert Chan Wrote:  lua> R = fn'r,s: s=r/(1+r); r*s - k/(r+s)'

I named s = (v/c) because I had used it previously, solving for r = -(1+q)

s = r/(1+r) = -(1+q)/(-q) = 1+1/q = (v/c)

We can rephrase R of r, in terms of s instead.
R guesses = cbrt(k/2) .. sqrt(k+1), can be carried over.

lua> S2 = fn's,r: r=s/(1-s); r*s - k/(r+s)'

S2 can accidentally shoot to "wrong" root, thus not recommended.
Unless with good guesses ... but S1 can use them too!