module Math.QuadEq where

import Control.Monad
import Math.Utils

{- |

\(ax^2+bx+c=0, a\neq 0\)

>>> solveQuadEq 1 (-3) 2
[1.0,2.0]
>>> solveQuadEq 1 (-2) 1
[1.0]
>>> solveQuadEq 1 0 0
[-0.0]
>>> solveQuadEq 1 0 (-2)
[-1.4142135623730951,1.414213562373095]
>>> solveQuadEq 1 (-1) (-1)
[-0.6180339887498948,1.618033988749895]
>>> solveQuadEq 1 (-2147483648) 2147483647
[1.0,2.147483647e9]
>>> solveQuadEq 1 (-2147483648) 1
[4.656612873077393e-10,2.147483648e9]
>>> solveQuadEq 0 1 1
*** Exception: solveQuadEq: (0,1,1)
-}
solveQuadEq :: Integer -> Integer -> Integer -> [Double]
solveQuadEq :: Integer -> Integer -> Integer -> [Double]
solveQuadEq Integer
a Integer
b Integer
c
  | Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = [Char] -> [Double]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [Double]) -> [Char] -> [Double]
forall a b. (a -> b) -> a -> b
$ [Char]
"solveQuadEq: " [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> (Integer, Integer, Integer) -> [Char]
forall a. Show a => a -> [Char]
show (Integer
a, Integer
b, Integer
c)
  | Integer
d Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = []
  | Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = [-Double
b' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a')]
  | Integer
b Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0
  , Double
x0 <- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (-Double
b' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt Double
d')
  , Double
x1 <- (-Double
b' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt Double
d') Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a') =
      if Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 then [Double
x0, Double
x1] else [Double
x1, Double
x0]
  | Double
x0 <- (-Double
b' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt Double
d') Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a')
  , Double
x1 <- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (-Double
b' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt Double
d') =
      if Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 then [Double
x0, Double
x1] else [Double
x1, Double
x0]
  where
    d :: Integer
d = Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
4 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
c
    a' :: Double
a' = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
a
    b' :: Double
b' = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
b
    c' :: Double
c' = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
c
    d' :: Double
d' = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
d

{- |

\(ax^2+bx+c=0, a\neq 0\)

>>> solveQuadEqInteger 1 (-3) 2
[1,2]
>>> solveQuadEqInteger 1 (-2) 1
[1]
>>> solveQuadEqInteger 1 0 0
[0]
>>> solveQuadEqInteger 1 0 (-2)
[]
>>> solveQuadEqInteger 1 (-1) (-1)
[]
>>> solveQuadEqInteger 0 1 1
*** Exception: solveQuadEqInteger: (0,1,1)
-}
solveQuadEqInteger :: Integer -> Integer -> Integer -> [Integer]
solveQuadEqInteger :: Integer -> Integer -> Integer -> [Integer]
solveQuadEqInteger Integer
a Integer
b Integer
c
  | Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = [Char] -> [Integer]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [Integer]) -> [Char] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [Char]
"solveQuadEqInteger: " [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> (Integer, Integer, Integer) -> [Char]
forall a. Show a => a -> [Char]
show (Integer
a, Integer
b, Integer
c)
  | Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = Integer -> Integer -> Integer -> [Integer]
solveQuadEqInteger (-Integer
a) (-Integer
b) (-Integer
c)
solveQuadEqInteger Integer
a Integer
b Integer
c
  | Integer
d Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = []
  | Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0, (Integer
x, Integer
0) <- Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem (-Integer
b) (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
a) = [Integer
x]
  | Bool
otherwise = do
      let !sqrtD :: Integer
sqrtD = Integer -> Integer
integerFloorSqrt Integer
d
      Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ Integer
sqrtD Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
sqrtD Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
d
      (!x, 0) <- (Integer -> (Integer, Integer))
-> [Integer] -> [(Integer, Integer)]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`quotRem` (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
a)) [-Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
sqrtD, -Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
sqrtD]
      pure x
  where
    !d :: Integer
d = Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
4 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
c