module Polynomial.Roots (roots) where
import Data.Complex
import Polynomial.Basic
roots :: RealFloat a => a
-> Int
-> [Complex a]
-> [Complex a]
roots :: forall a. RealFloat a => a -> Int -> [Complex a] -> [Complex a]
roots a
eps0 Int
count0 [Complex a]
as0 =
forall {t}.
RealFloat t =>
t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' a
eps0 Int
count0 [Complex a]
as0 []
where
roots' :: t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' t
eps Int
count [Complex t]
as [Complex t]
xs
| forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex t]
as forall a. Ord a => a -> a -> Bool
<= Int
2 = Complex t
xforall a. a -> [a] -> [a]
:[Complex t]
xs
| Bool
otherwise =
t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' t
eps Int
count (forall {a}. (Eq a, Num a) => a -> [a] -> [a] -> [a]
deflate Complex t
x [Complex t]
bs [forall a. [a] -> a
last [Complex t]
as]) (Complex t
xforall a. a -> [a] -> [a]
:[Complex t]
xs)
where
x :: Complex t
x = forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre t
eps Int
count [Complex t]
as Complex t
0
bs :: [Complex t]
bs = forall a. Int -> [a] -> [a]
drop Int
1 (forall a. [a] -> [a]
reverse (forall a. Int -> [a] -> [a]
drop Int
1 [Complex t]
as))
deflate :: a -> [a] -> [a] -> [a]
deflate a
z [a]
bs' [a]
cs
| [a]
bs' forall a. Eq a => a -> a -> Bool
== [] = [a]
cs
| Bool
otherwise =
a -> [a] -> [a] -> [a]
deflate a
z (forall a. [a] -> [a]
tail [a]
bs') (((forall a. [a] -> a
head [a]
bs')forall a. Num a => a -> a -> a
+a
zforall a. Num a => a -> a -> a
*(forall a. [a] -> a
head [a]
cs))forall a. a -> [a] -> [a]
:[a]
cs)
laguerre :: RealFloat a => a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre :: forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre a
eps0 Int
count [Complex a]
as0 Complex a
x0
| Int
count forall a. Ord a => a -> a -> Bool
<= Int
0 = Complex a
x0
| forall a. RealFloat a => Complex a -> a
magnitude (Complex a
x0 forall a. Num a => a -> a -> a
- Complex a
x0') forall a. Ord a => a -> a -> Bool
< a
eps0 = Complex a
x0'
| Bool
otherwise = forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre a
eps0 (Int
count forall a. Num a => a -> a -> a
- Int
1) [Complex a]
as0 Complex a
x0'
where x0' :: Complex a
x0' = forall {a}.
RealFloat a =>
a
-> [Complex a]
-> [Complex a]
-> [Complex a]
-> Complex a
-> Complex a
laguerre2 a
eps0 [Complex a]
as0 [Complex a]
as0' [Complex a]
as0'' Complex a
x0
as0' :: [Complex a]
as0' = forall a. Num a => [a] -> [a]
polyderiv [Complex a]
as0
as0'' :: [Complex a]
as0'' = forall a. Num a => [a] -> [a]
polyderiv [Complex a]
as0'
laguerre2 :: a
-> [Complex a]
-> [Complex a]
-> [Complex a]
-> Complex a
-> Complex a
laguerre2 a
eps [Complex a]
as [Complex a]
as' [Complex a]
as'' Complex a
x
| forall a. RealFloat a => Complex a -> a
magnitude Complex a
b forall a. Ord a => a -> a -> Bool
< a
eps = Complex a
x
| forall a. RealFloat a => Complex a -> a
magnitude Complex a
gp forall a. Ord a => a -> a -> Bool
< forall a. RealFloat a => Complex a -> a
magnitude Complex a
gm =
if Complex a
gm forall a. Eq a => a -> a -> Bool
== Complex a
0 then Complex a
x forall a. Num a => a -> a -> a
- Complex a
1 else Complex a
x forall a. Num a => a -> a -> a
- Complex a
nforall a. Fractional a => a -> a -> a
/Complex a
gm
| Bool
otherwise =
if Complex a
gp forall a. Eq a => a -> a -> Bool
== Complex a
0 then Complex a
x forall a. Num a => a -> a -> a
- Complex a
1 else Complex a
x forall a. Num a => a -> a -> a
- Complex a
nforall a. Fractional a => a -> a -> a
/Complex a
gp
where gp :: Complex a
gp = Complex a
g forall a. Num a => a -> a -> a
+ Complex a
delta
gm :: Complex a
gm = Complex a
g forall a. Num a => a -> a -> a
- Complex a
delta
g :: Complex a
g = Complex a
dforall a. Fractional a => a -> a -> a
/Complex a
b
delta :: Complex a
delta = forall a. Floating a => a -> a
sqrt ((Complex a
nforall a. Num a => a -> a -> a
-Complex a
1)forall a. Num a => a -> a -> a
*(Complex a
nforall a. Num a => a -> a -> a
*Complex a
h forall a. Num a => a -> a -> a
- Complex a
g2))
h :: Complex a
h = Complex a
g2 forall a. Num a => a -> a -> a
- Complex a
fforall a. Fractional a => a -> a -> a
/Complex a
b
b :: Complex a
b = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as Complex a
x
d :: Complex a
d = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as' Complex a
x
f :: Complex a
f = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as'' Complex a
x
g2 :: Complex a
g2 = Complex a
gforall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int)
n :: Complex a
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex a]
as)