module Diagrams.Solve.Polynomial
( quadForm
, cubForm
, quartForm
, cubForm'
, quartForm'
) where
import Data.List (maximumBy)
import Data.Ord (comparing)
import Prelude hiding ((^))
import qualified Prelude as P ((^))
tau :: Floating a => a
tau :: a
tau = 2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi
(^) :: (Num a) => a -> Integer -> a
^ :: a -> Integer -> a
(^) = a -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
(P.^)
aboutZero' :: (Ord a, Num a) => a -> a -> Bool
aboutZero' :: a -> a -> Bool
aboutZero' toler :: a
toler x :: a
x = a -> a
forall a. Num a => a -> a
abs a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
toler
{-# INLINE aboutZero' #-}
quadForm :: (Floating d, Ord d) => d -> d -> d -> [d]
quadForm :: d -> d -> d -> [d]
quadForm a :: d
a b :: d
b c :: d
c
| d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
&& d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
&& d
c d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = [0]
| d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
&& d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = []
| d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = [-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
b]
| d
d d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = []
| d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = [d -> d
forall a. Floating a => a -> a
sqrt (-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a), -d -> d
forall a. Floating a => a -> a
sqrt (-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a)]
| d
d d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = [-d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(2d -> d -> d
forall a. Num a => a -> a -> a
*d
a)]
| Bool
otherwise = [d
qd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a, d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
q]
where d :: d
d = d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2 d -> d -> d
forall a. Num a => a -> a -> a
- 4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
c
q :: d
q = -1d -> d -> d
forall a. Fractional a => a -> a -> a
/2d -> d -> d
forall a. Num a => a -> a -> a
*(d
b d -> d -> d
forall a. Num a => a -> a -> a
+ d -> d
forall a. Num a => a -> a
signum d
b d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
sqrt d
d)
{-# INLINE quadForm #-}
_quadForm_prop :: Double -> Double -> Double -> Bool
_quadForm_prop :: Double -> Double -> Double -> Bool
_quadForm_prop a :: Double
a b :: Double
b c :: Double
c = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' 1e-10 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm Double
a Double
b Double
c)
where eval :: Double -> Double
eval x :: Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c
cubForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
cubForm' :: d -> d -> d -> d -> d -> [d]
cubForm' toler :: d
toler a :: d
a b :: d
b c :: d
c d :: d
d
| d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
a = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm d
b d
c d
d
| d
delta d -> d -> Bool
forall a. Ord a => a -> a -> Bool
> 0 = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map d -> d
trig [0,1,2]
| d
delta d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
&& d
disc d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = [ -d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) ]
| d
delta d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
&& d
disc d -> d -> Bool
forall a. Eq a => a -> a -> Bool
/= 0 = [ (d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
- 9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
d)d -> d -> d
forall a. Fractional a => a -> a -> a
/(2d -> d -> d
forall a. Num a => a -> a -> a
*d
disc)
, (9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
- 4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
+ d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^3)d -> d -> d
forall a. Fractional a => a -> a -> a
/(d
a d -> d -> d
forall a. Num a => a -> a -> a
* d
disc)
]
| Bool
otherwise = [-d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) d -> d -> d
forall a. Num a => a -> a -> a
- d
ccd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) d -> d -> d
forall a. Num a => a -> a -> a
+ d
discd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cc)]
where delta :: d
delta = 18d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
- 4d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^3d -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
+ d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2 d -> d -> d
forall a. Num a => a -> a -> a
- 4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> Integer -> d
forall a. Num a => a -> Integer -> a
^3 d -> d -> d
forall a. Num a => a -> a -> a
- 27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Num a => a -> a -> a
*d
dd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2
disc :: d
disc = 3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
- d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2
qq :: d
qq = d -> d
forall a. Floating a => a -> a
sqrt(-27d -> d -> d
forall a. Num a => a -> a -> a
*(d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2)d -> d -> d
forall a. Num a => a -> a -> a
*d
delta)
qq' :: d
qq' = if d -> d
forall a. Num a => a -> a
abs (d
xx d -> d -> d
forall a. Num a => a -> a -> a
+ d
qq) d -> d -> Bool
forall a. Ord a => a -> a -> Bool
> d -> d
forall a. Num a => a -> a
abs (d
xx d -> d -> d
forall a. Num a => a -> a -> a
- d
qq) then d
qq else -d
qq
cc :: d
cc = d -> d
forall a. (Ord a, Floating a) => a -> a
cubert (1d -> d -> d
forall a. Fractional a => a -> a -> a
/2d -> d -> d
forall a. Num a => a -> a -> a
*(d
qq' d -> d -> d
forall a. Num a => a -> a -> a
+ d
xx))
xx :: d
xx = 2d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^3 d -> d -> d
forall a. Num a => a -> a -> a
- 9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
+ 27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Num a => a -> a -> a
*d
d
p :: d
p = d
discd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2)
q :: d
q = d
xxd -> d -> d
forall a. Fractional a => a -> a -> a
/(27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^3)
phi :: d
phi = 1d -> d -> d
forall a. Fractional a => a -> a -> a
/3d -> d -> d
forall a. Num a => a -> a -> a
*d -> d
forall a. Floating a => a -> a
acos(3d -> d -> d
forall a. Num a => a -> a -> a
*d
qd -> d -> d
forall a. Fractional a => a -> a -> a
/(2d -> d -> d
forall a. Num a => a -> a -> a
*d
p)d -> d -> d
forall a. Num a => a -> a -> a
*d -> d
forall a. Floating a => a -> a
sqrt(-3d -> d -> d
forall a. Fractional a => a -> a -> a
/d
p))
trig :: d -> d
trig k :: d
k = 2 d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
sqrt(-d
pd -> d -> d
forall a. Fractional a => a -> a -> a
/3) d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
cos(d
phi d -> d -> d
forall a. Num a => a -> a -> a
- d
kd -> d -> d
forall a. Num a => a -> a -> a
*d
forall a. Floating a => a
taud -> d -> d
forall a. Fractional a => a -> a -> a
/3) d -> d -> d
forall a. Num a => a -> a -> a
- d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(3d -> d -> d
forall a. Num a => a -> a -> a
*d
a)
cubert :: a -> a
cubert x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = -((-a
x)a -> a -> a
forall a. Floating a => a -> a -> a
**(1a -> a -> a
forall a. Fractional a => a -> a -> a
/3))
| Bool
otherwise = a
xa -> a -> a
forall a. Floating a => a -> a -> a
**(1a -> a -> a
forall a. Fractional a => a -> a -> a
/3)
{-# INLINE cubForm' #-}
cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm :: d -> d -> d -> d -> [d]
cubForm = d -> d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
cubForm' 1e-10
{-# INLINE cubForm #-}
_cubForm_prop :: Double -> Double -> Double -> Double -> Bool
_cubForm_prop :: Double -> Double -> Double -> Double -> Bool
_cubForm_prop a :: Double
a b :: Double
b c :: Double
c d :: Double
d = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' 1e-5 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm Double
a Double
b Double
c Double
d)
where eval :: Double -> Double
eval x :: Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d
quartForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d]
quartForm' :: d -> d -> d -> d -> d -> d -> [d]
quartForm' toler :: d
toler c4 :: d
c4 c3 :: d
c3 c2 :: d
c2 c1 :: d
c1 c0 :: d
c0
| d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
c4 = d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
c3 d
c2 d
c1 d
c0
| d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
c0 = 0 d -> [d] -> [d]
forall a. a -> [a] -> [a]
: d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
c4 d
c3 d
c2 d
c1
| Bool
otherwise = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map (\x :: d
x->d
xd -> d -> d
forall a. Num a => a -> a -> a
-(d
ad -> d -> d
forall a. Fractional a => a -> a -> a
/4)) [d]
roots
where
[a :: d
a,b :: d
b,c :: d
c,d :: d
d] = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map (d -> d -> d
forall a. Fractional a => a -> a -> a
/d
c4) [d
c3,d
c2,d
c1,d
c0]
p :: d
p = d
b d -> d -> d
forall a. Num a => a -> a -> a
- 3d -> d -> d
forall a. Fractional a => a -> a -> a
/8d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2
q :: d
q = 1d -> d -> d
forall a. Fractional a => a -> a -> a
/8d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^3d -> d -> d
forall a. Num a => a -> a -> a
-d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/2d -> d -> d
forall a. Num a => a -> a -> a
+d
c
r :: d
r = (-3d -> d -> d
forall a. Fractional a => a -> a -> a
/256)d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^4d -> d -> d
forall a. Num a => a -> a -> a
+d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/16d -> d -> d
forall a. Num a => a -> a -> a
-d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/4d -> d -> d
forall a. Num a => a -> a -> a
+d
d
roots :: [d]
roots | d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
r =
0 d -> [d] -> [d]
forall a. a -> [a] -> [a]
: d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm 1 0 d
p d
q
| d
u d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< -d
toler Bool -> Bool -> Bool
|| d
v d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< -d
toler = []
| Bool
otherwise = [d]
s1[d] -> [d] -> [d]
forall a. [a] -> [a] -> [a]
++[d]
s2
z :: d
z:_ = d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm 1 (-d
pd -> d -> d
forall a. Fractional a => a -> a -> a
/2) (-d
r) (d
pd -> d -> d
forall a. Num a => a -> a -> a
*d
rd -> d -> d
forall a. Fractional a => a -> a -> a
/2 d -> d -> d
forall a. Num a => a -> a -> a
- d
qd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2d -> d -> d
forall a. Fractional a => a -> a -> a
/8)
u :: d
u = d
zd -> Integer -> d
forall a. Num a => a -> Integer -> a
^2 d -> d -> d
forall a. Num a => a -> a -> a
- d
r
v :: d
v = 2d -> d -> d
forall a. Num a => a -> a -> a
*d
z d -> d -> d
forall a. Num a => a -> a -> a
- d
p
u' :: d
u' = if d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
u then 0 else d -> d
forall a. Floating a => a -> a
sqrt d
u
v' :: d
v' = if d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
v then 0 else d -> d
forall a. Floating a => a -> a
sqrt d
v
s1 :: [d]
s1 = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm 1 (if d
qd -> d -> Bool
forall a. Ord a => a -> a -> Bool
<0 then -d
v' else d
v') (d
zd -> d -> d
forall a. Num a => a -> a -> a
-d
u')
s2 :: [d]
s2 = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm 1 (if d
qd -> d -> Bool
forall a. Ord a => a -> a -> Bool
<0 then d
v' else -d
v') (d
zd -> d -> d
forall a. Num a => a -> a -> a
+d
u')
{-# INLINE quartForm' #-}
quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
quartForm :: d -> d -> d -> d -> d -> [d]
quartForm = d -> d -> d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d]
quartForm' 1e-10
{-# INLINE quartForm #-}
_quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool
_quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool
_quartForm_prop a :: Double
a b :: Double
b c :: Double
c d :: Double
d e :: Double
e = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' 1e-5 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
quartForm Double
a Double
b Double
c Double
d Double
e)
where eval :: Double -> Double
eval x :: Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e