Skip to main content
1 of 2
Eelvex
  • 5.5k
  • 1
  • 28
  • 43

#J, 222

247 - recursion + numerical

NB. multiple lines for visibility
f=:[:-(]$:~[:}.[-%&{.*],#&0@l)`[@.(0>l=:-&#)
r=:1={.
s=:[:=/+/@(2=&*/\])@(,.]*_1^i.@#)@({.@>@(([;];f;j;f f j=.]f f)([:}:(i._5)*])))@|.
u=:0,0([(]*_1^r@k)>:@]^:(1 1-:k=:[:*(0 p.~[)*(p.(,-)))^:_)~]
h=:{.@((n@],]{~[:r[:*p.*(p.n=:{.+-:@-~/))^:_ u)`('n'"_)@.s@|.

Methods

  • Using Sturm's theorem to find number of roots.

f is polynomial division (recursive)

s is Sturm's theorem output: 0 for no real roots, 1 otherwise. It evaluates Sturm's polynomials ([;];f;j;f f j=.]f f) and checks their behaviour at +/- infinity.

  • Using bisection method to avoid divisions by zero, local maxima etc.

u scans from 0 to +/- infinity to find an appropriate interval

((n@],]{~[:r[:*p.*(p.n=:{.+-:@-~/))^:_ u) evaluates the polynomial at the edges and midpoint then halfs the interval accordingly

Verification of results

Using J's buildin polynomial root finder as realRoots:

ti,>([;h;realRoots) each samples
┌──────────────┬────────┬─────────────────┐
│a b c d e     │output  │all real roots   │
├──────────────┼────────┼─────────────────┤
│_1 2 3 4 5    │_1.11029│3.37192 _1.11029 │
├──────────────┼────────┼─────────────────┤
│_1 2 3 4 _5   │0.728727│3.18248 0.728727 │
├──────────────┼────────┼─────────────────┤
│1 2 3 4 5     │n       │                 │
├──────────────┼────────┼─────────────────┤
│1 1 2 _1 _1   │0.728262│0.728262 _0.52236│
├──────────────┼────────┼─────────────────┤
│1 2 3 4 0     │0       │_1.65063 0       │
├──────────────┼────────┼─────────────────┤
│1 2 0 0 1     │_1      │_1.83929 _1      │
├──────────────┼────────┼─────────────────┤
│_1 _1 _1 _1 _1│n       │                 │
└──────────────┴────────┴─────────────────┘
Eelvex
  • 5.5k
  • 1
  • 28
  • 43