...the problem is that the variable r1 can become 0, which makes p2 and q2 empty (results of division by 0). Their absolute value is treated as 0, so the function accepts the wrong result as correct. This is why we don't see divergence at the tooltips, although there are a lot of them.
Here is the version of the function, which checks for this problems, too, and restarts the search. This seems to be robust enough if the coefficients remain below 1.e9. Just remove the tooltip lines if the code works as expected.
Code:
#SingleInstance Force
#NoEnv
SetBatchLines -1
SetFormat FloatFast, 0.17e ; maximum precision
CoordMode ToolTip, Screen
s = 1,0,0,0,0,
Random,,0
Loop {
Random r, -1.e6, 0
ToolTip %A_Index% %r%,180,0,1
PolyRoots(s r)
}
;------------------------------------------------------------------------------
PolyRoots(A) {
Loop Parse, A, `,
{
i := A_Index-1
If i = 0
an := A_loopfield
re%i% := im%i% := b%i% := d%i% := 0.0
a%A_Index% := A_loopfield / an
}
a0 := 0.0, b1 := d1 := 1.0, b2 := a2 ; for degree 1 poly
order := n := i, n1 := n+1, n2 := n+2
tol := 1.e-12, LIMIT := 1000
q1 := q2 := r1 := r2 := r3 := 0.0
i := j := k := cnt := 0 ; cnt = root counter
While n > 1 { ; The main Lin-Bairstow loop
p1 := q1 := .2 ; Initial values for factor p(x) = x^2 + p1*x + q1
Try:
Loop %LIMIT% {
i := 2, j := 1, k := 0
Loop %n% {
b%i% := a%i% - p1 * b%j% - q1 * b%k%
d%i% := b%i% - p1 * d%j% - q1 * d%k%
i++, j++, k++
}
j := n-1, k := n-2
r1 := d%j%**2 - (d%n% - b%n%) * d%k%
If abs(r1) < 1.e-99
break
p2 := (b%n% * d%j% - b%n1% * d%k%) / r1
q2 := (b%n1% * d%j% - (d%n% - b%n%) * b%n%) / r1
p1 := p1 + p2
q1 := q1 + q2
If (Abs(p2) < tol && Abs(q2) < tol)
GoTo Cont ; Converged!
}
ToolTip DIVERGED %A_Now%,0,0,2 ; For tests
Random p1, -0.5, 0.5
Random q1, -0.5, 0.5
GoTo Try ; If not converged, random re-try
Cont:
r1 := p1*p1 - 4*q1 ; Continue with recording the found roots
t := cnt + 1
If (r1 < 0) { ; imaginary roots
re%t% := re%cnt% :=-p1/2
im%t% :=-im%cnt% := sqrt(Abs(r1))/2
}
Else { ; real roots
r2 := sqrt(r1)
re%cnt% := (r2 - p1) / 2
re%t% := -(r2 + p1) / 2
im%t% :=im%cnt% := 0
}
cnt += 2 ; update root counter
n -= 2, n1 -= 2, n2 -= 2 ; reduce polynomial order
If (n > 1)
Loop %n2% ; coefficients of the new polynomial
a%A_Index% := b%A_Index%
}
if (n = 1) ; obtain last single real root
re%cnt% := -b2, im%cnt% := 0
Loop %order% ; format result for display in proportional font
i := A_index-1, j := round(im%i%,9), r := round(re%i%,9)
, x .= (r<0 ? "" : " ") r (j<0 ? " " : " +") j "i" (A_Index < order ? "`n" : "")
Return x
}
I also changed the format of the output for nicer look in MsgBox. You may want to replace multiple spaces with single ones, and "`n" with ",". There is no "+" in front of negative imaginary parts, either.
Login
Register
FAQ
Search
