{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE NoMonomorphismRestriction #-}

{- | Algebraic Real Numbers for exact computation

   Since 0.4.0.0
-}
module Algebra.Field.AlgebraicReal
  ( Algebraic,
    algebraic,

    -- * Operations
    nthRoot,
    nthRoot',
    improve,
    approximate,
    approxFractional,

    -- * Equation solver
    realRoots,
    complexRoots,

    -- * Interval arithmetic
    Interval (..),
    representative,
    includes,
    intersect,

    -- * Internal utility functions
    strum,
    presultant,
    factors,
    realPartPoly,
    imagPartPoly,
    sqFreePart,
  )
where

import Algebra.Prelude.Core hiding
  ( intersect,
    normalize,
  )
import Algebra.Ring.Polynomial.Factorise (clearDenom, factorHensel)
import Algebra.Ring.Polynomial.Univariate
import Control.Arrow ((***))
import Control.Lens
  ( each,
    ifoldMap,
    ix,
    (%~),
    (&),
    (*~),
    _Wrapped,
  )
import Control.Monad.Random
import Data.Bits
import qualified Data.Coerce as C
import qualified Data.IntMap as IM
import Data.Maybe (fromJust)
import Data.MonoTraversable
import qualified Data.Ratio as R
import qualified Data.Set as S
import qualified Data.Sized as SV
import GHC.Num (Num)
import Math.NumberTheory.Roots
import Numeric.Algebra.Complex hiding (i)
import Numeric.Decidable.Zero (isZero)
import Numeric.Domain.GCD (gcd)
import qualified Numeric.Field.Fraction as F
import Numeric.Semiring.ZeroProduct (ZeroProductSemiring)
import System.Entropy (getEntropy)
import System.IO.Unsafe (unsafePerformIO)
import qualified Prelude as P

stdGenFromEntropy :: IO StdGen
stdGenFromEntropy :: IO StdGen
stdGenFromEntropy = Int -> StdGen
mkStdGen (Int -> StdGen) -> (ByteString -> Int) -> ByteString -> StdGen
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Element ByteString -> Int -> Int) -> Int -> ByteString -> Int
forall mono b.
MonoFoldable mono =>
(Element mono -> b -> b) -> b -> mono -> b
ofoldr (\Element ByteString
a -> Int -> Int -> Int
forall r. Additive r => r -> r -> r
(+) (Word8 -> Int
forall a. Enum a => a -> Int
fromEnum Word8
Element ByteString
a) (Int -> Int) -> (Int -> Int) -> Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
8)) Int
0 (ByteString -> StdGen) -> IO ByteString -> IO StdGen
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> IO ByteString
getEntropy Int
8
{-# INLINE stdGenFromEntropy #-}

factors :: Unipol Rational -> [Unipol Rational]
factors :: Unipol Rational -> [Unipol Rational]
factors Unipol Rational
f =
  let (Integer
a, Unipol Integer
p) = Unipol Rational -> (Integer, Unipol Integer)
forall a.
(CoeffRing a, Euclidean a) =>
Unipol (Fraction a) -> (a, Unipol a)
clearDenom Unipol Rational
f
   in ((Int, Set (Unipol Integer)) -> [Unipol Rational])
-> [(Int, Set (Unipol Integer))] -> [Unipol Rational]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((Unipol Integer -> Unipol Rational)
-> [Unipol Integer] -> [Unipol Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Unipol Rational -> Unipol Rational
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize (Unipol Rational -> Unipol Rational)
-> (Unipol Integer -> Unipol Rational)
-> Unipol Integer
-> Unipol Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Rational) -> Unipol Integer -> Unipol Rational
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
a)) ([Unipol Integer] -> [Unipol Rational])
-> ((Int, Set (Unipol Integer)) -> [Unipol Integer])
-> (Int, Set (Unipol Integer))
-> [Unipol Rational]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Set (Unipol Integer) -> [Unipol Integer]
forall a. Set a -> [a]
S.toList (Set (Unipol Integer) -> [Unipol Integer])
-> ((Int, Set (Unipol Integer)) -> Set (Unipol Integer))
-> (Int, Set (Unipol Integer))
-> [Unipol Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Set (Unipol Integer)) -> Set (Unipol Integer)
forall a b. (a, b) -> b
snd) ([(Int, Set (Unipol Integer))] -> [Unipol Rational])
-> [(Int, Set (Unipol Integer))] -> [Unipol Rational]
forall a b. (a -> b) -> a -> b
$
        IntMap (Set (Unipol Integer)) -> [(Int, Set (Unipol Integer))]
forall a. IntMap a -> [(Int, a)]
IM.toList (IntMap (Set (Unipol Integer)) -> [(Int, Set (Unipol Integer))])
-> IntMap (Set (Unipol Integer)) -> [(Int, Set (Unipol Integer))]
forall a b. (a -> b) -> a -> b
$
          (Integer, IntMap (Set (Unipol Integer)))
-> IntMap (Set (Unipol Integer))
forall a b. (a, b) -> b
snd ((Integer, IntMap (Set (Unipol Integer)))
 -> IntMap (Set (Unipol Integer)))
-> (Integer, IntMap (Set (Unipol Integer)))
-> IntMap (Set (Unipol Integer))
forall a b. (a -> b) -> a -> b
$
            (Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
 -> StdGen -> (Integer, IntMap (Set (Unipol Integer))))
-> StdGen
-> Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
-> (Integer, IntMap (Set (Unipol Integer)))
forall a b c. (a -> b -> c) -> b -> a -> c
flip Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
-> StdGen -> (Integer, IntMap (Set (Unipol Integer)))
forall g a. Rand g a -> g -> a
evalRand (IO StdGen -> StdGen
forall a. IO a -> a
unsafePerformIO IO StdGen
stdGenFromEntropy) (Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
 -> (Integer, IntMap (Set (Unipol Integer))))
-> Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
-> (Integer, IntMap (Set (Unipol Integer)))
forall a b. (a -> b) -> a -> b
$
              Unipol Integer
-> Rand StdGen (Integer, IntMap (Set (Unipol Integer)))
forall (m :: * -> *).
MonadRandom m =>
Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
factorHensel Unipol Integer
p

{- | Algebraic real numbers, which can be expressed as a root
   of a rational polynomial.
-}
data Algebraic
  = Algebraic
      { Algebraic -> Unipol Rational
eqn :: Unipol Rational
      , Algebraic -> [Unipol Rational]
_strumseq :: [Unipol Rational]
      , Algebraic -> Interval Rational
_interval :: Interval Rational
      }
  | Rational !Rational

-- deriving (Show)
-- Invariants for constructor Algebraic:
--   (1) eqn must be monic and irreducible polynomial with degree > 1.
--   (2) strumseq is the standard strum sequence of eqn
--   (3) interval contains exactly one root of eqn

viewNthRoot :: Algebraic -> Maybe (Int, Rational)
viewNthRoot :: Algebraic -> Maybe (Int, Rational)
viewNthRoot (Algebraic Unipol Rational
eq [Unipol Rational]
_ Interval Rational
_) =
  let (Rational
lc, OrderedMonomial Grevlex 1
lm) = Unipol Rational
-> (Coefficient (Unipol Rational),
    OrderedMonomial
      (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
forall poly.
IsOrderedPolynomial poly =>
poly
-> (Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
leadingTerm Unipol Rational
eq
   in if Unipol Rational -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' (Unipol Rational
eq Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Group r => r -> r -> r
- (Coefficient (Unipol Rational),
 OrderedMonomial
   (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
-> Unipol Rational
forall poly.
IsOrderedPolynomial poly =>
(Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
-> poly
toPolynomial (Rational
Coefficient (Unipol Rational)
lc, OrderedMonomial
  (MOrder (Unipol Rational)) (Arity (Unipol Rational))
OrderedMonomial Grevlex 1
lm)) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then (Int, Rational) -> Maybe (Int, Rational)
forall a. a -> Maybe a
Just (OrderedMonomial Grevlex 1 -> Monomial 1
forall k (ordering :: k) (n :: Nat).
OrderedMonomial ordering n -> Monomial n
getMonomial OrderedMonomial Grevlex 1
lm Monomial 1 -> Ordinal 1 -> Int
forall (f :: * -> *) (n :: Nat) c.
(CFoldable f, Dom f c) =>
Sized f n c -> Ordinal n -> c
SV.%!! Ordinal 1
0, Rational -> Rational
forall r. Group r => r -> r
negate (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Coefficient (Unipol Rational)
forall poly. IsPolynomial poly => poly -> Coefficient poly
constantTerm Unipol Rational
eq)
        else Maybe (Int, Rational)
forall a. Maybe a
Nothing
viewNthRoot Algebraic
_ = Maybe (Int, Rational)
forall a. Maybe a
Nothing

showsNthRoot :: Int -> Rational -> ShowS
showsNthRoot :: Int -> Rational -> ShowS
showsNthRoot Int
1 Rational
q = Int -> Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
6 Rational
q
showsNthRoot Int
2 Rational
q = Char -> ShowS
showChar Char
'√' ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
6 Rational
q
showsNthRoot Int
n Rational
q = Int -> Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
6 Rational
q ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString String
"^(1/" ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> ShowS
forall a. Show a => a -> ShowS
shows Int
n ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
')'

instance Show Algebraic where
  showsPrec :: Int -> Algebraic -> ShowS
showsPrec Int
d (Rational Rational
q) = Int -> Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
d Rational
q
  showsPrec Int
d a :: Algebraic
a@(Algebraic Unipol Rational
eq [Unipol Rational]
str Interval Rational
int) = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
11) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
    case Algebraic -> Maybe (Int, Rational)
viewNthRoot Algebraic
a of
      Maybe (Int, Rational)
Nothing ->
        String -> ShowS
showString String
"Algebraic " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Unipol Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 Unipol Rational
eq ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
          ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Unipol Rational] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 [Unipol Rational]
str
          ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
          ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Interval Rational -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 Interval Rational
int
      Just (Int
nth, Rational
q)
        | Algebraic
a Algebraic -> Algebraic -> Bool
forall a. Ord a => a -> a -> Bool
< Algebraic
0 ->
          Bool -> ShowS -> ShowS
showParen (Int
5 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
d Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
11) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
            String -> ShowS
showString String
"-" ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Rational -> ShowS
showsNthRoot Int
nth Rational
q
        | Bool
otherwise -> Int -> Rational -> ShowS
showsNthRoot Int
nth Rational
q

negateA :: Algebraic -> Algebraic
negateA :: Algebraic -> Algebraic
negateA Algebraic
r = (-Integer
1 :: Integer) Integer -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
.* Algebraic
r
{-# INLINE [1] negateA #-}

minusA :: Algebraic -> Algebraic -> Algebraic
minusA :: Algebraic -> Algebraic -> Algebraic
minusA Algebraic
a Algebraic
b = Algebraic -> Algebraic -> Algebraic
plusA Algebraic
a (Algebraic -> Algebraic
negateA Algebraic
b)
{-# INLINE [1] minusA #-}

instance Additive Algebraic where
  + :: Algebraic -> Algebraic -> Algebraic
(+) = Algebraic -> Algebraic -> Algebraic
plusA
  {-# INLINE (+) #-}

instance Group Algebraic where
  negate :: Algebraic -> Algebraic
negate = Algebraic -> Algebraic
negateA
  {-# INLINE negate #-}
  (-) = Algebraic -> Algebraic -> Algebraic
minusA
  {-# INLINE (-) #-}

instance Monoidal Algebraic where
  zero :: Algebraic
zero = Rational -> Algebraic
Rational Rational
forall m. Monoidal m => m
zero
  {-# INLINE zero #-}

instance LeftModule Natural Algebraic where
  .* :: Natural -> Algebraic -> Algebraic
(.*) = Integer -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
(.*) (Integer -> Algebraic -> Algebraic)
-> (Natural -> Integer) -> Natural -> Algebraic -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Natural -> Integer
forall a. Integral a => a -> Integer
P.toInteger
  {-# INLINE (.*) #-}

instance RightModule Natural Algebraic where
  *. :: Algebraic -> Natural -> Algebraic
(*.) = (Natural -> Algebraic -> Algebraic)
-> Algebraic -> Natural -> Algebraic
forall a b c. (a -> b -> c) -> b -> a -> c
flip Natural -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
(.*)
  {-# INLINE (*.) #-}

instance LeftModule Integer Algebraic where
  .* :: Integer -> Algebraic -> Algebraic
(.*) = Algebraic -> Algebraic -> Algebraic
multA (Algebraic -> Algebraic -> Algebraic)
-> (Integer -> Algebraic) -> Integer -> Algebraic -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Algebraic
Rational (Rational -> Algebraic)
-> (Integer -> Rational) -> Integer -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral
  {-# INLINE (.*) #-}

instance LeftModule (Fraction Integer) Algebraic where
  .* :: Rational -> Algebraic -> Algebraic
(.*) = Algebraic -> Algebraic -> Algebraic
multA (Algebraic -> Algebraic -> Algebraic)
-> (Rational -> Algebraic) -> Rational -> Algebraic -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Algebraic
Rational
  {-# INLINE (.*) #-}

instance RightModule Integer Algebraic where
  *. :: Algebraic -> Integer -> Algebraic
(*.) = (Integer -> Algebraic -> Algebraic)
-> Algebraic -> Integer -> Algebraic
forall a b c. (a -> b -> c) -> b -> a -> c
flip Integer -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
(.*)
  {-# INLINE (*.) #-}

instance RightModule (Fraction Integer) Algebraic where
  *. :: Algebraic -> Rational -> Algebraic
(*.) = (Rational -> Algebraic -> Algebraic)
-> Algebraic -> Rational -> Algebraic
forall a b c. (a -> b -> c) -> b -> a -> c
flip Rational -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
(.*)
  {-# INLINE (*.) #-}

instance Multiplicative Algebraic where
  * :: Algebraic -> Algebraic -> Algebraic
(*) = Algebraic -> Algebraic -> Algebraic
multA
  {-# INLINE (*) #-}
  pow1p :: Algebraic -> Natural -> Algebraic
pow1p Algebraic
a = Algebraic -> Natural -> Algebraic
forall a. Integral a => Algebraic -> a -> Algebraic
powA Algebraic
a (Natural -> Algebraic)
-> (Natural -> Natural) -> Natural -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Natural -> Natural
forall a. Enum a => a -> a
succ
  {-# INLINE pow1p #-}

instance LeftModule (Scalar (Fraction Integer)) Algebraic where
  .* :: Scalar Rational -> Algebraic -> Algebraic
(.*) = Algebraic -> Algebraic -> Algebraic
multA (Algebraic -> Algebraic -> Algebraic)
-> (Scalar Rational -> Algebraic)
-> Scalar Rational
-> Algebraic
-> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Algebraic
Rational (Rational -> Algebraic)
-> (Scalar Rational -> Rational) -> Scalar Rational -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Scalar Rational -> Rational
forall r. Scalar r -> r
runScalar
  {-# INLINE (.*) #-}

instance RightModule (Scalar (Fraction Integer)) Algebraic where
  *. :: Algebraic -> Scalar Rational -> Algebraic
(*.) = (Scalar Rational -> Algebraic -> Algebraic)
-> Algebraic -> Scalar Rational -> Algebraic
forall a b c. (a -> b -> c) -> b -> a -> c
flip Scalar Rational -> Algebraic -> Algebraic
forall r m. LeftModule r m => r -> m -> m
(.*)
  {-# INLINE (*.) #-}

instance Eq Algebraic where
  Rational Rational
p == :: Algebraic -> Algebraic -> Bool
== Algebraic Unipol Rational
g [Unipol Rational]
_ Interval Rational
j =
    Rational -> Bool
forall r. DecidableZero r => r -> Bool
isZero (Scalar Rational -> Rational
forall r. Scalar r -> r
runScalar (Scalar Rational -> Rational) -> Scalar Rational -> Rational
forall a b. (a -> b) -> a -> b
$ (Ordinal (Arity (Unipol Rational)) -> Scalar Rational)
-> Unipol Rational -> Scalar Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Scalar Rational -> Ordinal 1 -> Scalar Rational
forall a b. a -> b -> a
const (Scalar Rational -> Ordinal 1 -> Scalar Rational)
-> Scalar Rational -> Ordinal 1 -> Scalar Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Scalar Rational
forall r. r -> Scalar r
Scalar Rational
p) Unipol Rational
g) Bool -> Bool -> Bool
&& Rational
p Rational -> Interval Rational -> Bool
forall a. Ord a => a -> Interval a -> Bool
`isIn` Interval Rational
j
  Algebraic Unipol Rational
g [Unipol Rational]
_ Interval Rational
j == Rational Rational
p =
    Rational -> Bool
forall r. DecidableZero r => r -> Bool
isZero (Scalar Rational -> Rational
forall r. Scalar r -> r
runScalar (Scalar Rational -> Rational) -> Scalar Rational -> Rational
forall a b. (a -> b) -> a -> b
$ (Ordinal (Arity (Unipol Rational)) -> Scalar Rational)
-> Unipol Rational -> Scalar Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Scalar Rational -> Ordinal 1 -> Scalar Rational
forall a b. a -> b -> a
const (Scalar Rational -> Ordinal 1 -> Scalar Rational)
-> Scalar Rational -> Ordinal 1 -> Scalar Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Scalar Rational
forall r. r -> Scalar r
Scalar Rational
p) Unipol Rational
g) Bool -> Bool -> Bool
&& Rational
p Rational -> Interval Rational -> Bool
forall a. Ord a => a -> Interval a -> Bool
`isIn` Interval Rational
j
  Rational Rational
p == Rational Rational
q = Rational
p Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
q
  Algebraic Unipol Rational
f [Unipol Rational]
ss Interval Rational
i == Algebraic Unipol Rational
g [Unipol Rational]
_ Interval Rational
j =
    Unipol Rational -> Unipol Rational
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize Unipol Rational
f Unipol Rational -> Unipol Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Unipol Rational -> Unipol Rational
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize Unipol Rational
g Bool -> Bool -> Bool
&& Interval Rational -> [Unipol Rational] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn (Interval Rational
i Interval Rational -> Interval Rational -> Interval Rational
forall a.
(Monoidal a, Ord a) =>
Interval a -> Interval a -> Interval a
`intersect` Interval Rational
j) [Unipol Rational]
ss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1

-- | Takes intersection of two intervals.
intersect :: (Monoidal a, Ord a) => Interval a -> Interval a -> Interval a
intersect :: Interval a -> Interval a -> Interval a
intersect Interval a
i Interval a
j
  | Interval a -> Interval a -> Bool
forall a. Ord a => Interval a -> Interval a -> Bool
disjoint Interval a
i Interval a
j = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval a
forall m. Monoidal m => m
zero a
forall m. Monoidal m => m
zero
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= Interval a -> a
forall r. Interval r -> r
lower Interval a
j = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
lower Interval a
j) (Interval a -> a
forall r. Interval r -> r
upper Interval a
i)
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> Interval a -> a
forall r. Interval r -> r
lower Interval a
i = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
lower Interval a
i) (Interval a -> a
forall r. Interval r -> r
upper Interval a
j)
intersect Interval a
_ Interval a
_ = String -> Interval a
forall a. HasCallStack => String -> a
error String
"intersect"

instance Ord Algebraic where
  compare :: Algebraic -> Algebraic -> Ordering
compare (Rational Rational
q) (Rational Rational
r) = Rational -> Rational -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Rational
q Rational
r
  compare (Rational Rational
q) (Algebraic Unipol Rational
f [Unipol Rational]
sf Interval Rational
i) =
    let i' :: Interval Rational
i' = (Interval Rational -> Bool)
-> (Interval Rational -> Interval Rational)
-> Interval Rational
-> Interval Rational
forall a. (a -> Bool) -> (a -> a) -> a -> a
until (Bool -> Bool
not (Bool -> Bool)
-> (Interval Rational -> Bool) -> Interval Rational -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational
q Rational -> Interval Rational -> Bool
forall a. Ord a => a -> Interval a -> Bool
`isIn`)) ([Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sf) Interval Rational
i
     in if Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Interval Rational -> Rational
forall r. Interval r -> r
lower Interval Rational
i'
          then Ordering
LT
          else
            if Unipol Rational -> Bool
forall r. DecidableZero r => r -> Bool
isZero (Unipol Rational
f Unipol Rational -> [Unipol Rational] -> Unipol Rational
forall poly (t :: * -> *).
(IsOrderedPolynomial poly, Field (Coefficient poly), Functor t,
 Foldable t) =>
poly -> t poly -> poly
`modPolynomial` [Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0 Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Group r => r -> r -> r
- Coefficient (Unipol Rational) -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff Rational
Coefficient (Unipol Rational)
q])
              then Ordering
EQ
              else Ordering
GT
  compare a :: Algebraic
a@(Algebraic Unipol Rational
_ [Unipol Rational]
_ Interval Rational
_) b :: Algebraic
b@(Rational Rational
_) = Ordering -> Ordering
flipOrd (Ordering -> Ordering) -> Ordering -> Ordering
forall a b. (a -> b) -> a -> b
$ Algebraic -> Algebraic -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Algebraic
b Algebraic
a
    where
      flipOrd :: Ordering -> Ordering
flipOrd Ordering
EQ = Ordering
EQ
      flipOrd Ordering
LT = Ordering
GT
      flipOrd Ordering
GT = Ordering
LT
  compare a :: Algebraic
a@(Algebraic Unipol Rational
_ [Unipol Rational]
sf Interval Rational
i) b :: Algebraic
b@(Algebraic Unipol Rational
_ [Unipol Rational]
sg Interval Rational
j)
    | Algebraic
a Algebraic -> Algebraic -> Bool
forall a. Eq a => a -> a -> Bool
== Algebraic
b = Ordering
EQ
    | Bool
otherwise =
      let (Interval Rational
i', Interval Rational
j') = ((Interval Rational, Interval Rational) -> Bool)
-> ((Interval Rational, Interval Rational)
    -> (Interval Rational, Interval Rational))
-> (Interval Rational, Interval Rational)
-> (Interval Rational, Interval Rational)
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((Interval Rational -> Interval Rational -> Bool)
-> (Interval Rational, Interval Rational) -> Bool
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Interval Rational -> Interval Rational -> Bool
forall a. Ord a => Interval a -> Interval a -> Bool
disjoint) ([Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sf (Interval Rational -> Interval Rational)
-> (Interval Rational -> Interval Rational)
-> (Interval Rational, Interval Rational)
-> (Interval Rational, Interval Rational)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** [Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sg) (Interval Rational
i, Interval Rational
j)
       in if Interval Rational -> Rational
forall r. Interval r -> r
upper Interval Rational
i' Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Interval Rational -> Rational
forall r. Interval r -> r
lower Interval Rational
j' then Ordering
LT else Ordering
GT

disjoint :: Ord a => Interval a -> Interval a -> Bool
disjoint :: Interval a -> Interval a -> Bool
disjoint Interval a
i Interval a
j = Interval a -> a
forall r. Interval r -> r
upper Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< Interval a -> a
forall r. Interval r -> r
lower Interval a
j Bool -> Bool -> Bool
|| Interval a -> a
forall r. Interval r -> r
upper Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< Interval a -> a
forall r. Interval r -> r
lower Interval a
i

instance Commutative Algebraic

instance Unital Algebraic where
  one :: Algebraic
one = Rational -> Algebraic
Rational Rational
1
  {-# INLINE one #-}

instance Semiring Algebraic

instance Abelian Algebraic

instance Rig Algebraic

instance Ring Algebraic

instance ZeroProductSemiring Algebraic

instance Division Algebraic where
  recip :: Algebraic -> Algebraic
recip = Algebraic -> Algebraic
recipA
  {-# INLINE recip #-}

instance DecidableZero Algebraic where
  isZero :: Algebraic -> Bool
isZero (Rational Rational
0) = Bool
True
  isZero Algebraic
_ = Bool
False
  {-# INLINE isZero #-}

instance DecidableUnits Algebraic where
  isUnit :: Algebraic -> Bool
isUnit (Rational Rational
0) = Bool
False
  isUnit Algebraic
_ = Bool
True
  {-# INLINE isUnit #-}

  recipUnit :: Algebraic -> Maybe Algebraic
recipUnit (Rational Rational
0) = Maybe Algebraic
forall a. Maybe a
Nothing
  recipUnit Algebraic
r = Algebraic -> Maybe Algebraic
forall a. a -> Maybe a
Just (Algebraic -> Maybe Algebraic) -> Algebraic -> Maybe Algebraic
forall a b. (a -> b) -> a -> b
$ Algebraic -> Algebraic
recipA Algebraic
r

instance DecidableAssociates Algebraic where
  isAssociate :: Algebraic -> Algebraic -> Bool
isAssociate (Rational Rational
0) (Rational Rational
0) = Bool
True
  isAssociate Algebraic {} Algebraic {} = Bool
True
  isAssociate Algebraic
_ Algebraic
_ = Bool
False
  {-# INLINE isAssociate #-}

instance UnitNormalForm Algebraic where
  splitUnit :: Algebraic -> (Algebraic, Algebraic)
splitUnit Algebraic
x =
    if Algebraic
x Algebraic -> Algebraic -> Bool
forall a. Eq a => a -> a -> Bool
== Algebraic
0
      then (Algebraic
0, Rational -> Algebraic
Rational Rational
1)
      else
        if Algebraic
x Algebraic -> Algebraic -> Bool
forall a. Ord a => a -> a -> Bool
< Algebraic
0
          then (Rational -> Algebraic
Rational (-Rational
1), Algebraic -> Algebraic
forall r. Group r => r -> r
negate Algebraic
x)
          else (Rational -> Algebraic
Rational Rational
1, Algebraic
x)
  {-# INLINE splitUnit #-}

instance P.Num Algebraic where
  + :: Algebraic -> Algebraic -> Algebraic
(+) = Algebraic -> Algebraic -> Algebraic
plusA
  {-# INLINE (+) #-}

  (-) = Algebraic -> Algebraic -> Algebraic
minusA
  {-# INLINE (-) #-}

  * :: Algebraic -> Algebraic -> Algebraic
(*) = Algebraic -> Algebraic -> Algebraic
multA
  {-# INLINE (*) #-}

  abs :: Algebraic -> Algebraic
abs = (WrapAlgebra Algebraic -> WrapAlgebra Algebraic)
-> Algebraic -> Algebraic
C.coerce (WrapAlgebra Algebraic -> WrapAlgebra Algebraic
forall a. Num a => a -> a
P.abs :: WrapAlgebra Algebraic -> WrapAlgebra Algebraic)
  {-# INLINE abs #-}

  signum :: Algebraic -> Algebraic
signum = (WrapAlgebra Algebraic -> WrapAlgebra Algebraic)
-> Algebraic -> Algebraic
C.coerce (WrapAlgebra Algebraic -> WrapAlgebra Algebraic
forall a. Num a => a -> a
P.signum :: WrapAlgebra Algebraic -> WrapAlgebra Algebraic)
  {-# INLINE signum #-}

  negate :: Algebraic -> Algebraic
negate = Algebraic -> Algebraic
negateA
  {-# INLINE negate #-}

  fromInteger :: Integer -> Algebraic
fromInteger = Rational -> Algebraic
Rational (Rational -> Algebraic)
-> (Integer -> Rational) -> Integer -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Rational
forall r. Num r => Integer -> r
fromInteger
  {-# INLINE fromInteger #-}

instance P.Fractional Algebraic where
  recip :: Algebraic -> Algebraic
recip = Algebraic -> Algebraic
recipA
  {-# INLINE recip #-}

  fromRational :: Rational -> Algebraic
fromRational Rational
r = Rational -> Algebraic
Rational (Rational -> Integer
forall a. Ratio a -> a
R.numerator Rational
r Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Rational -> Integer
forall a. Ratio a -> a
R.denominator Rational
r)
  {-# INLINE fromRational #-}

scale :: (Ord r, Multiplicative r, Monoidal r) => r -> Interval r -> Interval r
scale :: r -> Interval r -> Interval r
scale r
k (Interval r
l r
r)
  | r
k r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
forall m. Monoidal m => m
zero = r -> r -> Interval r
forall r. r -> r -> Interval r
Interval (r
k r -> r -> r
forall r. Multiplicative r => r -> r -> r
* r
r) (r
k r -> r -> r
forall r. Multiplicative r => r -> r -> r
* r
l)
  | Bool
otherwise = r -> r -> Interval r
forall r. r -> r -> Interval r
Interval (r
k r -> r -> r
forall r. Multiplicative r => r -> r -> r
* r
l) (r
k r -> r -> r
forall r. Multiplicative r => r -> r -> r
* r
r)

-- | Test if the former interval includes the latter.
includes :: Ord a => Interval a -> Interval a -> Bool
Interval a
i a
j includes :: Interval a -> Interval a -> Bool
`includes` Interval a
i' a
j' = a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
i' Bool -> Bool -> Bool
&& a
j' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
j
{-# INLINE includes #-}

instance Group r => Additive (Interval r) where
  + :: Interval r -> Interval r -> Interval r
(+) = Interval r -> Interval r -> Interval r
forall r. Group r => Interval r -> Interval r -> Interval r
plusInt
  {-# INLINE (+) #-}

plusInt :: Group r => Interval r -> Interval r -> Interval r
plusInt :: Interval r -> Interval r -> Interval r
plusInt (Interval r
l r
u) (Interval r
l' r
u') = r -> r -> Interval r
forall r. r -> r -> Interval r
Interval (r
l r -> r -> r
forall r. Additive r => r -> r -> r
+ r
l') (r
u r -> r -> r
forall r. Additive r => r -> r -> r
+ r
u')

rootSumPoly :: Unipol Rational -> Unipol Rational -> Unipol Rational
rootSumPoly :: Unipol Rational -> Unipol Rational -> Unipol Rational
rootSumPoly Unipol Rational
f Unipol Rational
g =
  Unipol Rational -> Unipol Rational
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
normalize (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
    Unipol (Unipol Rational)
-> Unipol (Unipol Rational) -> Unipol Rational
forall k. (Euclidean k, CoeffRing k) => Unipol k -> Unipol k -> k
presultant (Unipol Rational -> Unipol (Unipol Rational)
liftP Unipol Rational
f) (Unipol (Unipol Rational) -> Unipol Rational)
-> Unipol (Unipol Rational) -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
      (Ordinal (Arity (Unipol (Unipol Rational)))
 -> Unipol (Unipol Rational))
-> Unipol (Unipol Rational) -> Unipol (Unipol Rational)
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol (Unipol Rational) -> Ordinal 1 -> Unipol (Unipol Rational)
forall a b. a -> b -> a
const (Unipol (Unipol Rational) -> Ordinal 1 -> Unipol (Unipol Rational))
-> Unipol (Unipol Rational)
-> Ordinal 1
-> Unipol (Unipol Rational)
forall a b. (a -> b) -> a -> b
$ Coefficient (Unipol (Unipol Rational)) -> Unipol (Unipol Rational)
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff (Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0) Unipol (Unipol Rational)
-> Unipol (Unipol Rational) -> Unipol (Unipol Rational)
forall r. Group r => r -> r -> r
- Ordinal (Arity (Unipol (Unipol Rational)))
-> Unipol (Unipol Rational)
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol (Unipol Rational)))
0) (Unipol (Unipol Rational) -> Unipol (Unipol Rational))
-> Unipol (Unipol Rational) -> Unipol (Unipol Rational)
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Unipol (Unipol Rational)
liftP Unipol Rational
g
  where
    liftP :: Unipol Rational -> Unipol (Unipol Rational)
    liftP :: Unipol Rational -> Unipol (Unipol Rational)
liftP = (Rational -> Unipol Rational)
-> Unipol Rational -> Unipol (Unipol Rational)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol Rational -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff

sqFreePart ::
  (Eq r, Euclidean r, Division r) =>
  Unipol r ->
  Unipol r
sqFreePart :: Unipol r -> Unipol r
sqFreePart Unipol r
f = Unipol r
f Unipol r -> Unipol r -> Unipol r
forall d. Euclidean d => d -> d -> d
`quot` Unipol r -> Unipol r -> Unipol r
forall d. GCDDomain d => d -> d -> d
gcd Unipol r
f (Ordinal (Arity (Unipol r)) -> Unipol r -> Unipol r
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity (Unipol r))
0 Unipol r
f)

{-# RULES
"minus-cancel" [~1] forall x.
  minusA x x =
    Rational 0
"plusminus/right" [~1] forall x (y :: Algebraic).
  minusA (plusA x y) y =
    x
"plusminus/left" [~1] forall x (y :: Algebraic).
  minusA (plusA x y) x =
    y
"plusminus/left" [~1] forall x (y :: Algebraic).
  plusA x (minusA y x) =
    y
  #-}

plusA :: Algebraic -> Algebraic -> Algebraic
plusA :: Algebraic -> Algebraic -> Algebraic
plusA (Rational Rational
r) (Rational Rational
q) = Rational -> Algebraic
Rational (Rational
r Rational -> Rational -> Rational
forall r. Additive r => r -> r -> r
+ Rational
q)
plusA a :: Algebraic
a@(Algebraic Unipol Rational
_ [Unipol Rational]
_ Interval Rational
_) b :: Algebraic
b@(Rational Rational
_) = Algebraic -> Algebraic -> Algebraic
plusA Algebraic
b Algebraic
a
plusA (Rational Rational
r) (Algebraic Unipol Rational
f [Unipol Rational]
_ Interval Rational
i) =
  let f' :: Unipol Rational
f' = (Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0 Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Group r => r -> r -> r
- Coefficient (Unipol Rational) -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff Rational
Coefficient (Unipol Rational)
r) Unipol Rational
f
   in Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic) -> Maybe Algebraic -> Algebraic
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic Unipol Rational
f' (Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval (Interval Rational -> Rational
forall r. Interval r -> r
lower Interval Rational
i Rational -> Rational -> Rational
forall r. Additive r => r -> r -> r
+ Rational
r) (Interval Rational -> Rational
forall r. Interval r -> r
upper Interval Rational
i Rational -> Rational -> Rational
forall r. Additive r => r -> r -> r
+ Rational
r))
plusA a :: Algebraic
a@(Algebraic Unipol Rational
f [Unipol Rational]
_ Interval Rational
_) b :: Algebraic
b@(Algebraic Unipol Rational
g [Unipol Rational]
_ Interval Rational
_) =
  let fg :: Unipol Rational
fg = Unipol Rational -> Unipol Rational -> Unipol Rational
rootSumPoly Unipol Rational
f Unipol Rational
g
      iij :: Interval Rational
iij = (Interval Rational -> Interval Rational -> Interval Rational)
-> Unipol Rational -> Algebraic -> Algebraic -> Interval Rational
catcher Interval Rational -> Interval Rational -> Interval Rational
forall r. Group r => Interval r -> Interval r -> Interval r
plusInt Unipol Rational
fg Algebraic
a Algebraic
b
   in Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic) -> Maybe Algebraic -> Algebraic
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic Unipol Rational
fg Interval Rational
iij
{-# INLINE [1] plusA #-}

isIn :: Ord a => a -> Interval a -> Bool
a
x isIn :: a -> Interval a -> Bool
`isIn` Interval a
i a
j = a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
x Bool -> Bool -> Bool
&& a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
j

{- | Smart constructor.
   @'algebraic' f i@ represents the unique root of
   rational polynomial @f@ in the interval @i@.
   If no root is found, or more than one root belongs to
   the given interval, returns @'Nothing'@.
-}
algebraic :: Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic :: Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic Unipol Rational
f Interval Rational
i =
  let pss :: [(Int, Unipol Rational, [Unipol Rational])]
pss =
        [ (Int
n, Unipol Rational
p, [Unipol Rational]
ss)
        | Unipol Rational
p <- Unipol Rational -> [Unipol Rational]
factors (Unipol Rational -> [Unipol Rational])
-> Unipol Rational -> [Unipol Rational]
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Unipol Rational
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
sqFreePart Unipol Rational
f
        , let (Int
n, [Unipol Rational]
ss) = Unipol Rational -> Interval Rational -> (Int, [Unipol Rational])
countRootsIn Unipol Rational
p Interval Rational
i
        , Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
        ]
   in case [(Int, Unipol Rational, [Unipol Rational])]
pss of
        [(Int
1, Unipol Rational
p, [Unipol Rational]
ss)] ->
          Algebraic -> Maybe Algebraic
forall a. a -> Maybe a
Just (Algebraic -> Maybe Algebraic) -> Algebraic -> Maybe Algebraic
forall a b. (a -> b) -> a -> b
$
            if Unipol Rational -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Rational
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
              then Rational -> Algebraic
Rational (Rational -> Algebraic) -> Rational -> Algebraic
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall r. Group r => r -> r
negate (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Coefficient (Unipol Rational)
forall poly. IsPolynomial poly => poly -> Coefficient poly
constantTerm Unipol Rational
p
              else Unipol Rational
-> [Unipol Rational] -> Interval Rational -> Algebraic
Algebraic Unipol Rational
p [Unipol Rational]
ss Interval Rational
i
        [(Int, Unipol Rational, [Unipol Rational])]
_ -> Maybe Algebraic
forall a. Maybe a
Nothing

algebraic' :: Unipol Rational -> Interval Rational -> Algebraic
algebraic' :: Unipol Rational -> Interval Rational -> Algebraic
algebraic' Unipol Rational
p Interval Rational
i =
  if Unipol Rational -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Rational
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
    then Rational -> Algebraic
Rational (Rational -> Algebraic) -> Rational -> Algebraic
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall r. Group r => r -> r
negate (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Coefficient (Unipol Rational)
forall poly. IsPolynomial poly => poly -> Coefficient poly
constantTerm Unipol Rational
p
    else Unipol Rational
-> [Unipol Rational] -> Interval Rational -> Algebraic
Algebraic Unipol Rational
p (Unipol Rational -> [Unipol Rational]
strum Unipol Rational
p) Interval Rational
i

catcher ::
  (Interval Rational -> Interval Rational -> Interval Rational) ->
  Unipol Rational ->
  Algebraic ->
  Algebraic ->
  Interval Rational
catcher :: (Interval Rational -> Interval Rational -> Interval Rational)
-> Unipol Rational -> Algebraic -> Algebraic -> Interval Rational
catcher Interval Rational -> Interval Rational -> Interval Rational
app Unipol Rational
h (Algebraic Unipol Rational
_ [Unipol Rational]
sf Interval Rational
i0) (Algebraic Unipol Rational
_ [Unipol Rational]
sg Interval Rational
j0) =
  let lens :: [Rational]
lens = (Interval Rational -> Rational)
-> [Interval Rational] -> [Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map Interval Rational -> Rational
forall r. Group r => Interval r -> r
size ([Interval Rational] -> [Rational])
-> [Interval Rational] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> [Interval Rational]
isolateRoots Unipol Rational
h
      (Interval Rational
i', Interval Rational
j') =
        ((Interval Rational, Interval Rational) -> Bool)
-> ((Interval Rational, Interval Rational)
    -> (Interval Rational, Interval Rational))
-> (Interval Rational, Interval Rational)
-> (Interval Rational, Interval Rational)
forall a. (a -> Bool) -> (a -> a) -> a -> a
until
          (\(Interval Rational
i, Interval Rational
j) -> (Rational -> Bool) -> [Rational] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Interval Rational -> Rational
forall r. Group r => Interval r -> r
size (Interval Rational -> Interval Rational -> Interval Rational
app Interval Rational
i Interval Rational
j) Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<) [Rational]
lens)
          ([Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sf (Interval Rational -> Interval Rational)
-> (Interval Rational -> Interval Rational)
-> (Interval Rational, Interval Rational)
-> (Interval Rational, Interval Rational)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** [Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sg)
          (Interval Rational
i0, Interval Rational
j0)
   in Interval Rational
i' Interval Rational -> Interval Rational -> Interval Rational
`app` Interval Rational
j'
catcher Interval Rational -> Interval Rational -> Interval Rational
_ Unipol Rational
_ Algebraic
_ Algebraic
_ = String -> Interval Rational
forall a. HasCallStack => String -> a
error String
"rational is impossible"

normalize :: (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
normalize :: Unipol r -> Unipol r
normalize = Unipol r -> Unipol r
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize (Unipol r -> Unipol r)
-> (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol r -> Unipol r
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
sqFreePart

stabilize :: Algebraic -> Algebraic
stabilize :: Algebraic -> Algebraic
stabilize r :: Algebraic
r@Rational {} = Algebraic
r
stabilize ar :: Algebraic
ar@(Algebraic Unipol Rational
f [Unipol Rational]
ss Interval Rational
int)
  | Interval Rational -> Rational
forall r. Interval r -> r
lower Interval Rational
int Rational -> Rational -> Rational
forall r. Multiplicative r => r -> r -> r
* Interval Rational -> Rational
forall r. Interval r -> r
upper Interval Rational
int Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 = Algebraic
ar
  | Bool
otherwise =
    let lh :: Interval Rational
lh = Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval (Interval Rational -> Rational
forall r. Interval r -> r
lower Interval Rational
int) Rational
0
        uh :: Interval Rational
uh = Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval Rational
0 (Interval Rational -> Rational
forall r. Interval r -> r
upper Interval Rational
int)
     in Unipol Rational
-> [Unipol Rational] -> Interval Rational -> Algebraic
Algebraic Unipol Rational
f [Unipol Rational]
ss (Interval Rational -> Algebraic) -> Interval Rational -> Algebraic
forall a b. (a -> b) -> a -> b
$ if Interval Rational -> [Unipol Rational] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn Interval Rational
lh [Unipol Rational]
ss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Interval Rational
uh else Interval Rational
lh

-- | Unsafe version of @'nthRoot'@.
nthRoot' :: Int -> Algebraic -> Algebraic
nthRoot' :: Int -> Algebraic -> Algebraic
nthRoot' Int
n = Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic)
-> (Algebraic -> Maybe Algebraic) -> Algebraic -> Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Algebraic -> Maybe Algebraic
nthRoot Int
n

{- | @nthRoot n r@ tries to computes n-th root of
   the given algebraic real @r@.
   It returns @'Nothing'@ if it's undefined.

   See also @'nthRoot''@.
-}
nthRoot :: Int -> Algebraic -> Maybe Algebraic
nthRoot :: Int -> Algebraic -> Maybe Algebraic
nthRoot Int
1 Algebraic
a = Algebraic -> Maybe Algebraic
forall a. a -> Maybe a
Just Algebraic
a
nthRoot Int
n Algebraic
a
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = Algebraic -> Algebraic
recipA (Algebraic -> Algebraic) -> Maybe Algebraic -> Maybe Algebraic
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Algebraic -> Maybe Algebraic
nthRoot (Int -> Int
forall a. Num a => a -> a
abs Int
n) Algebraic
a
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Maybe Algebraic
forall a. Maybe a
Nothing
  | Int -> Bool
forall a. Integral a => a -> Bool
even Int
n Bool -> Bool -> Bool
&& Algebraic
a Algebraic -> Algebraic -> Bool
forall a. Ord a => a -> a -> Bool
< Algebraic
0 = Maybe Algebraic
forall a. Maybe a
Nothing
  | Bool
otherwise =
    case Algebraic
a of
      Rational Rational
p ->
        (Rational -> Maybe Algebraic)
-> (Interval Rational -> Maybe Algebraic)
-> Either Rational (Interval Rational)
-> Maybe Algebraic
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either
          (Algebraic -> Maybe Algebraic
forall a. a -> Maybe a
Just (Algebraic -> Maybe Algebraic)
-> (Rational -> Algebraic) -> Rational -> Maybe Algebraic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Algebraic
Rational)
          (Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic (Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0 Unipol Rational -> Natural -> Unipol Rational
forall r. Unital r => r -> Natural -> r
^ Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Group r => r -> r -> r
- Coefficient (Unipol Rational) -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff Rational
Coefficient (Unipol Rational)
p))
          (Either Rational (Interval Rational) -> Maybe Algebraic)
-> Either Rational (Interval Rational) -> Maybe Algebraic
forall a b. (a -> b) -> a -> b
$ Int -> Rational -> Either Rational (Interval Rational)
nthRootRat Int
n Rational
p
      Algebraic Unipol Rational
f [Unipol Rational]
_ Interval Rational
range ->
        Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic ((OrderedMonomial
   (MOrder (Unipol Rational)) (Arity (Unipol Rational))
 -> OrderedMonomial
      (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
-> Unipol Rational -> Unipol Rational
forall poly.
IsOrderedPolynomial poly =>
(OrderedMonomial (MOrder poly) (Arity poly)
 -> OrderedMonomial (MOrder poly) (Arity poly))
-> poly -> poly
mapMonomialMonotonic ((Monomial 1 -> Identity (Monomial 1))
-> OrderedMonomial Grevlex 1
-> Identity (OrderedMonomial Grevlex 1)
forall s t. Rewrapping s t => Iso s t (Unwrapped s) (Unwrapped t)
_Wrapped ((Monomial 1 -> Identity (Monomial 1))
 -> OrderedMonomial Grevlex 1
 -> Identity (OrderedMonomial Grevlex 1))
-> ((Int -> Identity Int) -> Monomial 1 -> Identity (Monomial 1))
-> (Int -> Identity Int)
-> OrderedMonomial Grevlex 1
-> Identity (OrderedMonomial Grevlex 1)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (Monomial 1)
-> Traversal' (Monomial 1) (IxValue (Monomial 1))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Index (Monomial 1)
0 ((Int -> Identity Int)
 -> OrderedMonomial Grevlex 1
 -> Identity (OrderedMonomial Grevlex 1))
-> Int -> OrderedMonomial Grevlex 1 -> OrderedMonomial Grevlex 1
forall a s t. Num a => ASetter s t a a -> a -> s -> t
*~ Int
2) Unipol Rational
f) (Int -> Interval Rational -> Interval Rational
nthRootInterval Int
n Interval Rational
range)

nthRootRat :: Int -> Rational -> Either Rational (Interval Rational)
nthRootRat :: Int -> Rational -> Either Rational (Interval Rational)
nthRootRat Int
1 Rational
r = Rational -> Either Rational (Interval Rational)
forall a b. a -> Either a b
Left Rational
r
nthRootRat Int
2 Rational
r =
  let (Integer
p, Integer
q) = (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r, Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
      (Integer
isp, Integer
isq) = (Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
p, Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
q)
   in case (Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactSquareRoot Integer
p, Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactSquareRoot Integer
q) of
        (Just Integer
p', Just Integer
q') -> Rational -> Either Rational (Interval Rational)
forall a b. a -> Either a b
Left (Integer
p' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q')
        (Maybe Integer
mp', Maybe Integer
mq') ->
          Interval Rational -> Either Rational (Interval Rational)
forall a b. b -> Either a b
Right (Interval Rational -> Either Rational (Interval Rational))
-> Interval Rational -> Either Rational (Interval Rational)
forall a b. (a -> b) -> a -> b
$
            Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isp Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isq Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mq')
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isp Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isq Maybe Integer
mq')
nthRootRat Int
3 Rational
r =
  let (Integer
p, Integer
q) = (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r, Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
      (Integer
isp, Integer
isq) = (Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot Integer
p, Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot Integer
q)
   in case (Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactCubeRoot Integer
p, Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactCubeRoot Integer
q) of
        (Just Integer
p', Just Integer
q') -> Rational -> Either Rational (Interval Rational)
forall a b. a -> Either a b
Left (Integer
p' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q')
        (Maybe Integer
mp', Maybe Integer
mq') ->
          Interval Rational -> Either Rational (Interval Rational)
forall a b. b -> Either a b
Right (Interval Rational -> Either Rational (Interval Rational))
-> Interval Rational -> Either Rational (Interval Rational)
forall a b. (a -> b) -> a -> b
$
            Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isp Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isq Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mq')
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isp Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isq Maybe Integer
mq')
nthRootRat Int
4 Rational
r =
  let (Integer
p, Integer
q) = (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r, Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
      (Integer
isp, Integer
isq) = (Integer -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Integer
4 Integer
p, Integer -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Integer
4 Integer
q)
   in case (Integer -> Integer -> Maybe Integer
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Integer
4 Integer
p, Integer -> Integer -> Maybe Integer
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Integer
4 Integer
q) of
        (Just Integer
p', Just Integer
q') -> Rational -> Either Rational (Interval Rational)
forall a b. a -> Either a b
Left (Integer
p' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q')
        (Maybe Integer
mp', Maybe Integer
mq') ->
          Interval Rational -> Either Rational (Interval Rational)
forall a b. b -> Either a b
Right (Interval Rational -> Either Rational (Interval Rational))
-> Interval Rational -> Either Rational (Interval Rational)
forall a b. (a -> b) -> a -> b
$
            Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isp Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isq Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mq')
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isp Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isq Maybe Integer
mq')
nthRootRat Int
n Rational
r =
  let (Integer
p, Integer
q) = (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r, Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
      (Integer
isp, Integer
isq) = (Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n Integer
p, Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n Integer
q)
   in case (Int -> Integer -> Maybe Integer
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Int
n Integer
p, Int -> Integer -> Maybe Integer
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Int
n Integer
q) of
        (Just Integer
p', Just Integer
q') -> Rational -> Either Rational (Interval Rational)
forall a b. a -> Either a b
Left (Integer
p' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q')
        (Maybe Integer
mp', Maybe Integer
mq') ->
          Interval Rational -> Either Rational (Interval Rational)
forall a b. b -> Either a b
Right (Interval Rational -> Either Rational (Interval Rational))
-> Interval Rational -> Either Rational (Interval Rational)
forall a b. (a -> b) -> a -> b
$
            Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isp Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isq Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mq')
              (Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe (Integer
isp Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1) Maybe Integer
mp' Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer -> Maybe Integer -> Integer
forall a. a -> Maybe a -> a
fromMaybe Integer
isq Maybe Integer
mq')

nthRootRatCeil :: Int -> Rational -> Rational
nthRootRatCeil :: Int -> Rational -> Rational
nthRootRatCeil Int
1 Rational
r = Rational
r
nthRootRatCeil Int
2 Rational
r =
  let p :: Integer
p = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
      q :: Integer
q = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatCeil Int
3 Rational
r =
  let p :: Integer
p = Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
      q :: Integer
q = Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatCeil Int
4 Rational
r =
  let p :: Integer
p = Integer -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Integer
4 (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
      q :: Integer
q = Integer -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Integer
4 (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatCeil Int
n Rational
r =
  let p :: Integer
p = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
      q :: Integer
q = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r)
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q

nthRootRatFloor :: Int -> Rational -> Rational
nthRootRatFloor :: Int -> Rational -> Rational
nthRootRatFloor Int
1 Rational
r = Rational
r
nthRootRatFloor Int
2 Rational
r =
  let p :: Integer
p = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r)
      q :: Integer
q = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatFloor Int
3 Rational
r =
  let p :: Integer
p = Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r)
      q :: Integer
q = Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatFloor Int
4 Rational
r =
  let p :: Integer
p = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot (Int
4 :: Int) (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r)
      q :: Integer
q = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot (Int
4 :: Int) (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q
nthRootRatFloor Int
n Rational
r =
  let p :: Integer
p = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n (Rational -> Integer
forall t. Fraction t -> t
F.numerator Rational
r)
      q :: Integer
q = Int -> Integer -> Integer
forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot Int
n (Rational -> Integer
forall t. Fraction t -> t
F.denominator Rational
r) Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
1
   in Integer
p Integer -> Integer -> Rational
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
q

nthRootInterval :: Int -> Interval Rational -> Interval Rational
nthRootInterval :: Int -> Interval Rational -> Interval Rational
nthRootInterval Int
0 Interval Rational
_ = String -> Interval Rational
forall a. HasCallStack => String -> a
error String
"0-th root????????"
nthRootInterval Int
1 Interval Rational
i = Interval Rational
i
nthRootInterval Int
n (Interval Rational
lb Rational
ub)
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = Interval Rational -> Interval Rational
forall r. (Num r, Ord r, Division r) => Interval r -> Interval r
recipInt (Interval Rational -> Interval Rational)
-> Interval Rational -> Interval Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval (Int -> Rational -> Rational
nthRootRatFloor (Int -> Int
forall a. Num a => a -> a
abs Int
n) Rational
lb) (Int -> Rational -> Rational
nthRootRatCeil (Int -> Int
forall a. Num a => a -> a
abs Int
n) Rational
ub)
  | Bool
otherwise = Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval (Int -> Rational -> Rational
nthRootRatFloor Int
n Rational
lb) (Int -> Rational -> Rational
nthRootRatCeil Int
n Rational
ub)

rootMultPoly :: Unipol Rational -> Unipol Rational -> Unipol Rational
rootMultPoly :: Unipol Rational -> Unipol Rational -> Unipol Rational
rootMultPoly Unipol Rational
f Unipol Rational
g =
  let ts :: Map
  (OrderedMonomial
     (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
  (Coefficient (Unipol Rational))
ts = Unipol Rational
-> Map
     (OrderedMonomial
        (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
     (Coefficient (Unipol Rational))
forall poly.
IsOrderedPolynomial poly =>
poly
-> Map
     (OrderedMonomial (MOrder poly) (Arity poly)) (Coefficient poly)
terms Unipol Rational
g
      d :: Int
d = Unipol Rational -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Rational
g
      g' :: Unipol (Unipol Rational)
g' = Add (Unipol (Unipol Rational)) -> Unipol (Unipol Rational)
forall a. Add a -> a
runAdd (Add (Unipol (Unipol Rational)) -> Unipol (Unipol Rational))
-> Add (Unipol (Unipol Rational)) -> Unipol (Unipol Rational)
forall a b. (a -> b) -> a -> b
$ (OrderedMonomial Grevlex 1
 -> Rational -> Add (Unipol (Unipol Rational)))
-> Map (OrderedMonomial Grevlex 1) Rational
-> Add (Unipol (Unipol Rational))
forall i (f :: * -> *) m a.
(FoldableWithIndex i f, Monoid m) =>
(i -> a -> m) -> f a -> m
ifoldMap OrderedMonomial Grevlex 1
-> Rational -> Add (Unipol (Unipol Rational))
loop Map
  (OrderedMonomial
     (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
  (Coefficient (Unipol Rational))
Map (OrderedMonomial Grevlex 1) Rational
ts
      loop :: OrderedMonomial Grevlex 1
-> Rational -> Add (Unipol (Unipol Rational))
loop OrderedMonomial Grevlex 1
m Rational
k =
        let deg :: Int
deg = OrderedMonomial Grevlex 1 -> Int
forall k (ord :: k) (n :: Nat). OrderedMonomial ord n -> Int
totalDegree OrderedMonomial Grevlex 1
m
         in Unipol (Unipol Rational) -> Add (Unipol (Unipol Rational))
forall a. a -> Add a
Add (Unipol (Unipol Rational) -> Add (Unipol (Unipol Rational)))
-> Unipol (Unipol Rational) -> Add (Unipol (Unipol Rational))
forall a b. (a -> b) -> a -> b
$
              (Coefficient (Unipol (Unipol Rational)),
 OrderedMonomial
   (MOrder (Unipol (Unipol Rational)))
   (Arity (Unipol (Unipol Rational))))
-> Unipol (Unipol Rational)
forall poly.
IsOrderedPolynomial poly =>
(Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
-> poly
toPolynomial
                ( Coefficient (Unipol Rational) -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff Rational
Coefficient (Unipol Rational)
k Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Multiplicative r => r -> r -> r
* IsLabel "x" (Unipol Rational)
Unipol Rational
#x Unipol Rational -> Natural -> Unipol Rational
forall r. Unital r => r -> Natural -> r
^ Int -> Natural
forall a b. (Integral a, Num b) => a -> b
P.fromIntegral Int
deg
                , Monomial 1 -> OrderedMonomial Grevlex 1
forall k (ordering :: k) (n :: Nat).
Monomial n -> OrderedMonomial ordering n
OrderedMonomial (Monomial 1 -> OrderedMonomial Grevlex 1)
-> Monomial 1 -> OrderedMonomial Grevlex 1
forall a b. (a -> b) -> a -> b
$ Int -> Monomial 1
forall (f :: * -> *) a. (CPointed f, Dom f a) => a -> Sized f 1 a
singleton (Int -> Monomial 1) -> Int -> Monomial 1
forall a b. (a -> b) -> a -> b
$ Int
d Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
deg
                )
   in Unipol (Unipol Rational)
-> Unipol (Unipol Rational) -> Unipol Rational
forall k. (Euclidean k, CoeffRing k) => Unipol k -> Unipol k -> k
presultant (Unipol Rational -> Unipol (Unipol Rational)
liftP Unipol Rational
f) Unipol (Unipol Rational)
g'
  where
    liftP :: Unipol Rational -> Unipol (Unipol Rational)
    liftP :: Unipol Rational -> Unipol (Unipol Rational)
liftP = (Rational -> Unipol Rational)
-> Unipol Rational -> Unipol (Unipol Rational)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol Rational -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff

instance (Ord r, Multiplicative r, Monoidal r) => Multiplicative (Interval r) where
  * :: Interval r -> Interval r -> Interval r
(*) = Interval r -> Interval r -> Interval r
forall a.
(Monoidal a, Ord a, Multiplicative a) =>
Interval a -> Interval a -> Interval a
multInt
  {-# INLINE (*) #-}

multInt :: (Monoidal a, Ord a, Multiplicative a) => Interval a -> Interval a -> Interval a
multInt :: Interval a -> Interval a -> Interval a
multInt Interval a
i Interval a
j
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
forall m. Monoidal m => m
zero Bool -> Bool -> Bool
&& Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
forall m. Monoidal m => m
zero = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
upper Interval a
i a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
upper Interval a
j) (Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
lower Interval a
j)
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
forall m. Monoidal m => m
zero Bool -> Bool -> Bool
&& Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
forall m. Monoidal m => m
zero = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
upper Interval a
i a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
lower Interval a
j) (Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
upper Interval a
j)
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
forall m. Monoidal m => m
zero Bool -> Bool -> Bool
&& Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
forall m. Monoidal m => m
zero = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
upper Interval a
j a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
lower Interval a
i) (Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
upper Interval a
i)
  | Interval a -> a
forall r. Interval r -> r
lower Interval a
i a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
forall m. Monoidal m => m
zero Bool -> Bool -> Bool
&& Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
forall m. Monoidal m => m
zero = a -> a -> Interval a
forall r. r -> r -> Interval r
Interval (Interval a -> a
forall r. Interval r -> r
lower Interval a
j a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
lower Interval a
i) (Interval a -> a
forall r. Interval r -> r
upper Interval a
j a -> a -> a
forall r. Multiplicative r => r -> r -> r
* Interval a -> a
forall r. Interval r -> r
upper Interval a
i)
multInt Interval a
_ Interval a
_ = String -> Interval a
forall a. HasCallStack => String -> a
error String
"multInt"
{-# INLINEABLE multInt #-}

multA :: Algebraic -> Algebraic -> Algebraic
multA :: Algebraic -> Algebraic -> Algebraic
multA (Rational Rational
0) Algebraic
_ = Rational -> Algebraic
Rational Rational
0
multA Algebraic
_ (Rational Rational
0) = Rational -> Algebraic
Rational Rational
0
multA (Rational Rational
a) (Rational Rational
b) = Rational -> Algebraic
Rational (Rational
a Rational -> Rational -> Rational
forall r. Multiplicative r => r -> r -> r
* Rational
b)
multA a :: Algebraic
a@Rational {} b :: Algebraic
b@Algebraic {} = Algebraic -> Algebraic -> Algebraic
multA Algebraic
b Algebraic
a
multA (Algebraic Unipol Rational
f [Unipol Rational]
_ Interval Rational
i) (Rational Rational
a) =
  let factor :: Unipol Rational -> Unipol Rational
factor = Unipol Rational -> Unipol Rational
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize (Unipol Rational -> Unipol Rational)
-> (Unipol Rational -> Unipol Rational)
-> Unipol Rational
-> Unipol Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall r. Division r => r -> r
recip Rational
a Rational -> Unipol Rational -> Unipol Rational
forall r m. Module (Scalar r) m => r -> m -> m
.*. Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0)
   in Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic) -> Maybe Algebraic -> Algebraic
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic (Unipol Rational -> Unipol Rational
factor Unipol Rational
f) (Rational -> Interval Rational -> Interval Rational
forall r.
(Ord r, Multiplicative r, Monoidal r) =>
r -> Interval r -> Interval r
scale Rational
a Interval Rational
i)
-- ? Can't we use the strum sequence given beforehand?
multA Algebraic
a Algebraic
b =
  let fg :: Unipol Rational
fg = Unipol Rational -> Unipol Rational -> Unipol Rational
rootMultPoly (Algebraic -> Unipol Rational
eqn Algebraic
a) (Algebraic -> Unipol Rational
eqn Algebraic
b)
      int :: Interval Rational
int = (Interval Rational -> Interval Rational -> Interval Rational)
-> Unipol Rational -> Algebraic -> Algebraic -> Interval Rational
catcher Interval Rational -> Interval Rational -> Interval Rational
forall a.
(Monoidal a, Ord a, Multiplicative a) =>
Interval a -> Interval a -> Interval a
multInt Unipol Rational
fg (Algebraic -> Algebraic
stabilize Algebraic
a) (Algebraic -> Algebraic
stabilize Algebraic
b)
   in Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic) -> Maybe Algebraic -> Algebraic
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic Unipol Rational
fg Interval Rational
int

improveNonzero :: Algebraic -> Interval Rational
improveNonzero :: Algebraic -> Interval Rational
improveNonzero (Algebraic Unipol Rational
_ [Unipol Rational]
ss Interval Rational
int0) = Interval Rational -> Interval Rational
go Interval Rational
int0
  where
    go :: Interval Rational -> Interval Rational
go Interval Rational
int =
      let (Interval Rational
lh, Interval Rational
bh) = Interval Rational -> (Interval Rational, Interval Rational)
forall r.
(Ring r, Division r) =>
Interval r -> (Interval r, Interval r)
bisect Interval Rational
int
       in if Interval Rational -> [Unipol Rational] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn Interval Rational
lh [Unipol Rational]
ss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
            then Interval Rational
bh
            else Interval Rational -> Interval Rational
go Interval Rational
lh
improveNonzero Algebraic
_ = String -> Interval Rational
forall a. HasCallStack => String -> a
error String
"improveNonzero: Rational"

recipInt :: (Num r, Ord r, Division r) => Interval r -> Interval r
recipInt :: Interval r -> Interval r
recipInt Interval r
i
  | Interval r -> r
forall r. Interval r -> r
lower Interval r
i r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
0 = r -> r -> Interval r
forall r. r -> r -> Interval r
Interval (r -> r
forall r. Division r => r -> r
recip (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ Interval r -> r
forall r. Interval r -> r
upper Interval r
i) (r -> r
forall r. Division r => r -> r
recip (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ Interval r -> r
forall r. Interval r -> r
lower Interval r
i)
  | Bool
otherwise = r -> r -> Interval r
forall r. r -> r -> Interval r
Interval (r -> r
forall r. Division r => r -> r
recip (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ Interval r -> r
forall r. Interval r -> r
lower Interval r
i) (r -> r
forall r. Division r => r -> r
recip (r -> r) -> r -> r
forall a b. (a -> b) -> a -> b
$ Interval r -> r
forall r. Interval r -> r
upper Interval r
i)

recipA :: Algebraic -> Algebraic
recipA :: Algebraic -> Algebraic
recipA (Rational Rational
a) = Rational -> Algebraic
Rational (Rational -> Rational
forall r. Division r => r -> r
recip Rational
a)
recipA a :: Algebraic
a@Algebraic {} =
  let fi :: Unipol Rational
fi = Unipol Rational -> Unipol Rational
rootRecipPoly (Algebraic -> Unipol Rational
eqn Algebraic
a)
      sf :: [Unipol Rational]
sf = Unipol Rational -> [Unipol Rational]
strum Unipol Rational
fi
      i0 :: Interval Rational
i0 = Algebraic -> Interval Rational
improveNonzero Algebraic
a
      ls :: [Rational]
ls = (Interval Rational -> Rational)
-> [Interval Rational] -> [Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map Interval Rational -> Rational
forall r. Group r => Interval r -> r
size ([Interval Rational] -> [Rational])
-> [Interval Rational] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> [Interval Rational]
isolateRoots Unipol Rational
fi
      i' :: Interval Rational
i' = (Interval Rational -> Bool)
-> (Interval Rational -> Interval Rational)
-> Interval Rational
-> Interval Rational
forall a. (a -> Bool) -> (a -> a) -> a -> a
until (\Interval Rational
i -> (Rational -> Bool) -> [Rational] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Interval Rational -> Rational
forall r. Group r => Interval r -> r
size (Interval Rational -> Interval Rational
forall r. (Num r, Ord r, Division r) => Interval r -> Interval r
recipInt Interval Rational
i)) [Rational]
ls) ([Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
sf) Interval Rational
i0
   in Maybe Algebraic -> Algebraic
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Algebraic -> Algebraic) -> Maybe Algebraic -> Algebraic
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Interval Rational -> Maybe Algebraic
algebraic Unipol Rational
fi (Interval Rational -> Maybe Algebraic)
-> Interval Rational -> Maybe Algebraic
forall a b. (a -> b) -> a -> b
$ Interval Rational -> Interval Rational
forall r. (Num r, Ord r, Division r) => Interval r -> Interval r
recipInt Interval Rational
i'

rootRecipPoly :: Unipol Rational -> Unipol Rational
rootRecipPoly :: Unipol Rational -> Unipol Rational
rootRecipPoly Unipol Rational
f =
  let ts :: Map
  (OrderedMonomial
     (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
  (Coefficient (Unipol Rational))
ts = Unipol Rational
-> Map
     (OrderedMonomial
        (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
     (Coefficient (Unipol Rational))
forall poly.
IsOrderedPolynomial poly =>
poly
-> Map
     (OrderedMonomial (MOrder poly) (Arity poly)) (Coefficient poly)
terms Unipol Rational
f
      d :: Int
d = Unipol Rational -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Rational
f
   in Add (Unipol Rational) -> Unipol Rational
forall a. Add a -> a
runAdd (Add (Unipol Rational) -> Unipol Rational)
-> Add (Unipol Rational) -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
        (OrderedMonomial Grevlex 1 -> Rational -> Add (Unipol Rational))
-> Map (OrderedMonomial Grevlex 1) Rational
-> Add (Unipol Rational)
forall i (f :: * -> *) m a.
(FoldableWithIndex i f, Monoid m) =>
(i -> a -> m) -> f a -> m
ifoldMap
          ( \OrderedMonomial Grevlex 1
m Rational
k ->
              Unipol Rational -> Add (Unipol Rational)
forall a. a -> Add a
Add (Unipol Rational -> Add (Unipol Rational))
-> Unipol Rational -> Add (Unipol Rational)
forall a b. (a -> b) -> a -> b
$
                (Coefficient (Unipol Rational),
 OrderedMonomial
   (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
-> Unipol Rational
forall poly.
IsOrderedPolynomial poly =>
(Coefficient poly, OrderedMonomial (MOrder poly) (Arity poly))
-> poly
toPolynomial
                  (Rational
Coefficient (Unipol Rational)
k, Monomial 1 -> OrderedMonomial Grevlex 1
forall k (ordering :: k) (n :: Nat).
Monomial n -> OrderedMonomial ordering n
OrderedMonomial (Monomial 1 -> OrderedMonomial Grevlex 1)
-> Monomial 1 -> OrderedMonomial Grevlex 1
forall a b. (a -> b) -> a -> b
$ Int -> Monomial 1
forall (f :: * -> *) a. (CPointed f, Dom f a) => a -> Sized f 1 a
singleton (Int -> Monomial 1) -> Int -> Monomial 1
forall a b. (a -> b) -> a -> b
$ Int
d Int -> Int -> Int
forall r. Group r => r -> r -> r
- OrderedMonomial Grevlex 1 -> Int
forall k (ord :: k) (n :: Nat). OrderedMonomial ord n -> Int
totalDegree OrderedMonomial Grevlex 1
m)
          )
          Map
  (OrderedMonomial
     (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
  (Coefficient (Unipol Rational))
Map (OrderedMonomial Grevlex 1) Rational
ts

strum :: Unipol Rational -> [Unipol Rational]
strum :: Unipol Rational -> [Unipol Rational]
strum Unipol Rational
f =
  (Unipol Rational -> Unipol Rational -> Unipol Rational)
-> [Unipol Rational] -> [Unipol Rational] -> [Unipol Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Multiplicative r => r -> r -> r
(*) ([Unipol Rational] -> [Unipol Rational]
forall a. [a] -> [a]
cycle [Unipol Rational
1, Unipol Rational
1, -Unipol Rational
1, -Unipol Rational
1]) ([Unipol Rational] -> [Unipol Rational])
-> [Unipol Rational] -> [Unipol Rational]
forall a b. (a -> b) -> a -> b
$
    ((Unipol Rational, Unipol Rational, Unipol Rational)
 -> Unipol Rational)
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
-> [Unipol Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\(Unipol Rational
p, Unipol Rational
_, Unipol Rational
_) -> Unipol Rational
p Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Multiplicative r => r -> r -> r
* Coefficient (Unipol Rational) -> Unipol Rational
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff (Rational -> Rational
forall r. Division r => r -> r
recip (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall a. Num a => a -> a
abs (Unipol Rational -> Coefficient (Unipol Rational)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Rational
p))) ([(Unipol Rational, Unipol Rational, Unipol Rational)]
 -> [Unipol Rational])
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
-> [Unipol Rational]
forall a b. (a -> b) -> a -> b
$
      [(Unipol Rational, Unipol Rational, Unipol Rational)]
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
forall a. [a] -> [a]
reverse ([(Unipol Rational, Unipol Rational, Unipol Rational)]
 -> [(Unipol Rational, Unipol Rational, Unipol Rational)])
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
forall a b. (a -> b) -> a -> b
$ Unipol Rational
-> Unipol Rational
-> [(Unipol Rational, Unipol Rational, Unipol Rational)]
forall r. Euclidean r => r -> r -> [(r, r, r)]
prs Unipol Rational
f (Ordinal (Arity (Unipol Rational))
-> Unipol Rational -> Unipol Rational
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity (Unipol Rational))
0 Unipol Rational
f)

data Interval r = Interval {Interval r -> r
lower :: !r, Interval r -> r
upper :: !r} deriving (Interval r -> Interval r -> Bool
(Interval r -> Interval r -> Bool)
-> (Interval r -> Interval r -> Bool) -> Eq (Interval r)
forall r. Eq r => Interval r -> Interval r -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Interval r -> Interval r -> Bool
$c/= :: forall r. Eq r => Interval r -> Interval r -> Bool
== :: Interval r -> Interval r -> Bool
$c== :: forall r. Eq r => Interval r -> Interval r -> Bool
Eq, Eq (Interval r)
Eq (Interval r)
-> (Interval r -> Interval r -> Ordering)
-> (Interval r -> Interval r -> Bool)
-> (Interval r -> Interval r -> Bool)
-> (Interval r -> Interval r -> Bool)
-> (Interval r -> Interval r -> Bool)
-> (Interval r -> Interval r -> Interval r)
-> (Interval r -> Interval r -> Interval r)
-> Ord (Interval r)
Interval r -> Interval r -> Bool
Interval r -> Interval r -> Ordering
Interval r -> Interval r -> Interval r
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall r. Ord r => Eq (Interval r)
forall a. Ord a => Interval a -> Interval a -> Bool
forall r. Ord r => Interval r -> Interval r -> Ordering
forall r. Ord r => Interval r -> Interval r -> Interval r
min :: Interval r -> Interval r -> Interval r
$cmin :: forall r. Ord r => Interval r -> Interval r -> Interval r
max :: Interval r -> Interval r -> Interval r
$cmax :: forall r. Ord r => Interval r -> Interval r -> Interval r
>= :: Interval r -> Interval r -> Bool
$c>= :: forall a. Ord a => Interval a -> Interval a -> Bool
> :: Interval r -> Interval r -> Bool
$c> :: forall a. Ord a => Interval a -> Interval a -> Bool
<= :: Interval r -> Interval r -> Bool
$c<= :: forall a. Ord a => Interval a -> Interval a -> Bool
< :: Interval r -> Interval r -> Bool
$c< :: forall a. Ord a => Interval a -> Interval a -> Bool
compare :: Interval r -> Interval r -> Ordering
$ccompare :: forall r. Ord r => Interval r -> Interval r -> Ordering
$cp1Ord :: forall r. Ord r => Eq (Interval r)
Ord)

size :: Group r => Interval r -> r
size :: Interval r -> r
size (Interval r
l r
r) = r
r r -> r -> r
forall r. Group r => r -> r -> r
- r
l

instance Show r => Show (Interval r) where
  showsPrec :: Int -> Interval r -> ShowS
showsPrec Int
_ (Interval r
l r
u) = Char -> ShowS
showChar Char
'(' ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> r -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
5 r
l ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString String
", " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> r -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
5 r
u ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
']'

bisect :: (Ring r, Division r) => Interval r -> (Interval r, Interval r)
bisect :: Interval r -> (Interval r, Interval r)
bisect (Interval r
l r
u) =
  let mid :: r
mid = (r
l r -> r -> r
forall r. Additive r => r -> r -> r
+ r
u) r -> r -> r
forall r. Division r => r -> r -> r
/ Integer -> r
forall r. Ring r => Integer -> r
fromInteger' Integer
2
   in (r -> r -> Interval r
forall r. r -> r -> Interval r
Interval r
l r
mid, r -> r -> Interval r
forall r. r -> r -> Interval r
Interval r
mid r
u)

signChange :: (Ord a, Ring a, DecidableZero a) => [a] -> Int
signChange :: [a] -> Int
signChange [a]
xs =
  let nzs :: [a]
nzs = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> (a -> Bool) -> a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Bool
forall r. DecidableZero r => r -> Bool
isZero) [a]
xs
   in [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> Int) -> [a] -> Int
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
forall m. Monoidal m => m
zero) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall r. Multiplicative r => r -> r -> r
(*) [a]
nzs ([a] -> [a]
forall a. [a] -> [a]
tail [a]
nzs)

countRootsIn :: Unipol Rational -> Interval Rational -> (Int, [Unipol Rational])
countRootsIn :: Unipol Rational -> Interval Rational -> (Int, [Unipol Rational])
countRootsIn Unipol Rational
f Interval Rational
ints =
  let ss :: [Unipol Rational]
ss = Unipol Rational -> [Unipol Rational]
strum Unipol Rational
f
   in (Interval Rational -> [Unipol Rational] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn Interval Rational
ints [Unipol Rational]
ss, [Unipol Rational]
ss)

countChangeIn ::
  (Ord a, CoeffRing a) =>
  Interval a ->
  [Unipol a] ->
  Int
countChangeIn :: Interval a -> [Unipol a] -> Int
countChangeIn Interval a
ints [Unipol a]
ss =
  a -> [Unipol a] -> Int
forall a. (Ord a, CoeffRing a) => a -> [Unipol a] -> Int
signChangeAt (Interval a -> a
forall r. Interval r -> r
lower Interval a
ints) [Unipol a]
ss Int -> Int -> Int
forall r. Group r => r -> r -> r
- a -> [Unipol a] -> Int
forall a. (Ord a, CoeffRing a) => a -> [Unipol a] -> Int
signChangeAt (Interval a -> a
forall r. Interval r -> r
upper Interval a
ints) [Unipol a]
ss

signChangeAt :: (Ord a, CoeffRing a) => a -> [Unipol a] -> Int
signChangeAt :: a -> [Unipol a] -> Int
signChangeAt a
a [Unipol a]
fs = [a] -> Int
forall a. (Ord a, Ring a, DecidableZero a) => [a] -> Int
signChange ([a] -> Int) -> [a] -> Int
forall a b. (a -> b) -> a -> b
$ (Unipol a -> a) -> [Unipol a] -> [a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Scalar a -> a
forall r. Scalar r -> r
runScalar (Scalar a -> a) -> (Unipol a -> Scalar a) -> Unipol a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ordinal (Arity (Unipol a)) -> Scalar a) -> Unipol a -> Scalar a
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Scalar a -> Ordinal 1 -> Scalar a
forall a b. a -> b -> a
const (Scalar a -> Ordinal 1 -> Scalar a)
-> Scalar a -> Ordinal 1 -> Scalar a
forall a b. (a -> b) -> a -> b
$ a -> Scalar a
forall r. r -> Scalar r
Scalar a
a)) [Unipol a]
fs

rootBound ::
  (Num r, Ord r, CoeffRing r, Division r, UnitNormalForm r) =>
  Unipol r ->
  r
rootBound :: Unipol r -> r
rootBound Unipol r
f
  | Unipol r -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol r
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = r
0
  | Bool
otherwise =
    let a :: Coefficient (Unipol r)
a = Unipol r -> Coefficient (Unipol r)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol r
f
     in r
1 r -> r -> r
forall r. Additive r => r -> r -> r
+ Map (Monomial 1) r -> r
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((r -> r) -> Map (Monomial 1) r -> Map (Monomial 1) r
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (r -> r
forall r. UnitNormalForm r => r -> r
normaliseUnit (r -> r) -> (r -> r) -> r -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (r -> r -> r
forall r. Division r => r -> r -> r
/ r
Coefficient (Unipol r)
a)) (Map (Monomial 1) r -> Map (Monomial 1) r)
-> Map (Monomial 1) r -> Map (Monomial 1) r
forall a b. (a -> b) -> a -> b
$ Unipol r
-> Map (Monomial (Arity (Unipol r))) (Coefficient (Unipol r))
forall poly.
IsPolynomial poly =>
poly -> Map (Monomial (Arity poly)) (Coefficient poly)
terms' Unipol r
f)

isolateRoots :: Unipol Rational -> [Interval Rational]
isolateRoots :: Unipol Rational -> [Interval Rational]
isolateRoots Unipol Rational
f =
  let bd :: Rational
bd = Unipol Rational -> Rational
forall r.
(Num r, Ord r, CoeffRing r, Division r, UnitNormalForm r) =>
Unipol r -> r
rootBound Unipol Rational
f Rational -> Rational -> Rational
forall r. Multiplicative r => r -> r -> r
* Rational
2
   in Interval Rational -> [Interval Rational]
go (Rational -> Rational -> Interval Rational
forall r. r -> r -> Interval r
Interval (- Rational
bd) Rational
bd)
  where
    !ss :: [Unipol Rational]
ss = Unipol Rational -> [Unipol Rational]
strum Unipol Rational
f
    go :: Interval Rational -> [Interval Rational]
go Interval Rational
int =
      let rcount :: Int
rcount = Interval Rational -> [Unipol Rational] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn Interval Rational
int [Unipol Rational]
ss
       in if Int
rcount Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
            then []
            else
              if Int
rcount Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
                then [Interval Rational
int]
                else
                  let (Interval Rational
ls, Interval Rational
us) = Interval Rational -> (Interval Rational, Interval Rational)
forall r.
(Ring r, Division r) =>
Interval r -> (Interval r, Interval r)
bisect Interval Rational
int
                   in Interval Rational -> [Interval Rational]
go Interval Rational
ls [Interval Rational] -> [Interval Rational] -> [Interval Rational]
forall w. Monoid w => w -> w -> w
++ Interval Rational -> [Interval Rational]
go Interval Rational
us

{- | @'improve' r@ returns the same algebraic number,
   but with more tighter bounds.
-}
improve :: Algebraic -> Algebraic
improve :: Algebraic -> Algebraic
improve (Algebraic Unipol Rational
f [Unipol Rational]
ss Interval Rational
int) = Unipol Rational
-> [Unipol Rational] -> Interval Rational -> Algebraic
Algebraic Unipol Rational
f [Unipol Rational]
ss (Interval Rational -> Algebraic) -> Interval Rational -> Algebraic
forall a b. (a -> b) -> a -> b
$ [Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
ss Interval Rational
int
improve Algebraic
a = Algebraic
a

improveWith ::
  (Ord a, CoeffRing a, Division a) =>
  [Unipol a] ->
  Interval a ->
  Interval a
improveWith :: [Unipol a] -> Interval a -> Interval a
improveWith [Unipol a]
ss Interval a
int =
  let (Interval a
ls, Interval a
us) = Interval a -> (Interval a, Interval a)
forall r.
(Ring r, Division r) =>
Interval r -> (Interval r, Interval r)
bisect Interval a
int
   in if Interval a -> [Unipol a] -> Int
forall a. (Ord a, CoeffRing a) => Interval a -> [Unipol a] -> Int
countChangeIn Interval a
ls [Unipol a]
ss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Interval a
us else Interval a
ls

powA :: Integral a => Algebraic -> a -> Algebraic
powA :: Algebraic -> a -> Algebraic
powA Algebraic
r a
n
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = Int -> Algebraic -> Algebraic
nthRoot' (Int -> Int
forall a. Num a => a -> a
abs (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n) Algebraic
r
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = Rational -> Algebraic
Rational Rational
1
  | Bool
otherwise = Algebraic
r Algebraic -> Natural -> Algebraic
forall r. Unital r => r -> Natural -> r
^ a -> Natural
forall a b. (Integral a, Num b) => a -> b
P.fromIntegral a
n
{-# INLINE powA #-}

iterateImproveStrum :: Rational -> [Unipol Rational] -> Interval Rational -> Interval Rational
iterateImproveStrum :: Rational
-> [Unipol Rational] -> Interval Rational -> Interval Rational
iterateImproveStrum Rational
eps [Unipol Rational]
ss Interval Rational
int0 =
  (Interval Rational -> Bool)
-> (Interval Rational -> Interval Rational)
-> Interval Rational
-> Interval Rational
forall a. (a -> Bool) -> (a -> a) -> a -> a
until (\Interval Rational
int -> Interval Rational -> Rational
forall r. Group r => Interval r -> r
size Interval Rational
int Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
eps) ([Unipol Rational] -> Interval Rational -> Interval Rational
forall a.
(Ord a, CoeffRing a, Division a) =>
[Unipol a] -> Interval a -> Interval a
improveWith [Unipol Rational]
ss) Interval Rational
int0

fromFraction :: P.Fractional a => Fraction Integer -> a
fromFraction :: Rational -> a
fromFraction Rational
q = Integer -> a
forall r. Num r => Integer -> r
P.fromInteger (Rational -> Integer
forall t. Fraction t -> t
numerator Rational
q) a -> a -> a
forall a. Fractional a => a -> a -> a
P./ Integer -> a
forall r. Num r => Integer -> r
P.fromInteger (Rational -> Integer
forall t. Fraction t -> t
denominator Rational
q)

-- | Pseudo resultant. should we expose this?
presultant ::
  (Euclidean k, CoeffRing k) =>
  Unipol k ->
  Unipol k ->
  k
presultant :: Unipol k -> Unipol k -> k
presultant = Coefficient (Unipol k)
-> Coefficient (Unipol k)
-> Unipol k
-> Unipol k
-> Coefficient (Unipol k)
forall poly.
(Euclidean (Coefficient poly), IsOrderedPolynomial poly) =>
Coefficient poly
-> Coefficient poly -> poly -> poly -> Coefficient poly
go Coefficient (Unipol k)
forall r. Unital r => r
one Coefficient (Unipol k)
forall r. Unital r => r
one
  where
    go :: Coefficient poly
-> Coefficient poly -> poly -> poly -> Coefficient poly
go !Coefficient poly
res !Coefficient poly
acc poly
h poly
s
      | poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 =
        let (poly
_, poly
r) = poly
h poly -> poly -> (poly, poly)
forall k poly.
(k ~ Coefficient poly, Euclidean k, IsOrderedPolynomial poly) =>
poly -> poly -> (poly, poly)
`pDivModPoly` poly
s
            (Int
l, Int
k, Int
m) = (poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
h, poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
r, poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
s)
            res' :: Coefficient poly
res' =
              Coefficient poly
res Coefficient poly -> Coefficient poly -> Coefficient poly
forall r. Multiplicative r => r -> r -> r
* Int -> Coefficient poly
forall r n. (Ring r, Integral n) => n -> r
powSign (Int
l Int -> Int -> Int
forall r. Multiplicative r => r -> r -> r
* Int
m)
                Coefficient poly -> Coefficient poly -> Coefficient poly
forall r. Multiplicative r => r -> r -> r
* Coefficient poly -> Natural -> Coefficient poly
forall r. Unital r => r -> Natural -> r
pow (poly -> Coefficient poly
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff poly
s) (Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Int -> Int
forall a. Num a => a -> a
abs (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
l Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
k)
            adj :: Coefficient poly
adj = Coefficient poly -> Natural -> Coefficient poly
forall r. Unital r => r -> Natural -> r
pow (poly -> Coefficient poly
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff poly
s) (Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Int -> Int
forall a. Num a => a -> a
abs (Int
1 Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
l Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
m))
         in Coefficient poly
-> Coefficient poly -> poly -> poly -> Coefficient poly
go Coefficient poly
res' (Coefficient poly
acc Coefficient poly -> Coefficient poly -> Coefficient poly
forall r. Multiplicative r => r -> r -> r
* Coefficient poly
adj) poly
s poly
r
      | poly -> Bool
forall r. DecidableZero r => r -> Bool
isZero poly
h Bool -> Bool -> Bool
|| poly -> Bool
forall r. DecidableZero r => r -> Bool
isZero poly
s = Coefficient poly
forall m. Monoidal m => m
zero
      | poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
h Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = Coefficient poly -> Natural -> Coefficient poly
forall r. Unital r => r -> Natural -> r
pow (poly -> Coefficient poly
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff poly
s) (Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ poly -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' poly
h) Coefficient poly -> Coefficient poly -> Coefficient poly
forall r. Multiplicative r => r -> r -> r
* Coefficient poly
res Coefficient poly -> Coefficient poly -> Coefficient poly
forall d. Euclidean d => d -> d -> d
`quot` Coefficient poly
acc
      | Bool
otherwise = Coefficient poly
res Coefficient poly -> Coefficient poly -> Coefficient poly
forall d. Euclidean d => d -> d -> d
`quot` Coefficient poly
acc

powSign :: (Ring r, Integral n) => n -> r
powSign :: n -> r
powSign n
n
  | n -> Bool
forall a. Integral a => a -> Bool
even n
n = r
forall r. Unital r => r
one
  | Bool
otherwise = r -> r
forall r. Group r => r -> r
negate r
forall r. Unital r => r
one
{-# INLINE powSign #-}

-- | @'realRoots' f@ finds all real roots of the rational polynomial @f@.
realRoots :: Unipol Rational -> [Algebraic]
realRoots :: Unipol Rational -> [Algebraic]
realRoots = Unipol Rational -> [Algebraic]
realRootsIrreducible (Unipol Rational -> [Algebraic])
-> (Unipol Rational -> Unipol Rational)
-> Unipol Rational
-> [Algebraic]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol Rational -> Unipol Rational
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize (Unipol Rational -> [Algebraic])
-> (Unipol Rational -> [Unipol Rational])
-> Unipol Rational
-> [Algebraic]
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< Unipol Rational -> [Unipol Rational]
factors (Unipol Rational -> [Unipol Rational])
-> (Unipol Rational -> Unipol Rational)
-> Unipol Rational
-> [Unipol Rational]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol Rational -> Unipol Rational
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
sqFreePart

{- | Same as @'realRoots'@, but assumes that the given polynomial
   is monic and irreducible.
-}
realRootsIrreducible :: Unipol Rational -> [Algebraic]
realRootsIrreducible :: Unipol Rational -> [Algebraic]
realRootsIrreducible Unipol Rational
f =
  [ Unipol Rational -> Interval Rational -> Algebraic
algebraic' Unipol Rational
f Interval Rational
i
  | Interval Rational
i <- Unipol Rational -> [Interval Rational]
isolateRoots Unipol Rational
f
  ]

instance InvolutiveMultiplication Algebraic where
  adjoint :: Algebraic -> Algebraic
adjoint = Algebraic -> Algebraic
forall a. a -> a
id

instance TriviallyInvolutive Algebraic

{- | @'realRoots' f@ finds all complex roots of the rational polynomial @f@.

 CAUTION: This function currently comes with really naive implementation.
 Easy to explode.
-}
complexRoots :: Unipol Rational -> [Complex Algebraic]
complexRoots :: Unipol Rational -> [Complex Algebraic]
complexRoots = Unipol Rational -> [Complex Algebraic]
complexRoots' (Unipol Rational -> [Complex Algebraic])
-> (Unipol Rational -> [Unipol Rational])
-> Unipol Rational
-> [Complex Algebraic]
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< Unipol Rational -> [Unipol Rational]
factors (Unipol Rational -> [Unipol Rational])
-> (Unipol Rational -> Unipol Rational)
-> Unipol Rational
-> [Unipol Rational]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol Rational -> Unipol Rational
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
sqFreePart

complexRoots' :: Unipol Rational -> [Complex Algebraic]
complexRoots' :: Unipol Rational -> [Complex Algebraic]
complexRoots' Unipol Rational
f =
  let rp :: Unipol Rational
rp = Unipol Rational -> Unipol Rational
realPartPoly Unipol Rational
f
      ip :: Unipol Rational
ip = Unipol Rational -> Unipol Rational
imagPartPoly Unipol Rational
f
   in [ Complex Algebraic
c
      | Algebraic
r <- Unipol Rational -> [Algebraic]
realRoots Unipol Rational
rp
      , Algebraic
i <- Unipol Rational -> [Algebraic]
realRoots Unipol Rational
ip
      , let c :: Complex Algebraic
c = Algebraic -> Algebraic -> Complex Algebraic
forall a. a -> a -> Complex a
Complex Algebraic
r Algebraic
i
      , (Ordinal (Arity (Unipol Rational)) -> Complex Algebraic)
-> Unipol Rational -> Complex Algebraic
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Complex Algebraic -> Ordinal 1 -> Complex Algebraic
forall a b. a -> b -> a
const (Complex Algebraic -> Ordinal 1 -> Complex Algebraic)
-> Complex Algebraic -> Ordinal 1 -> Complex Algebraic
forall a b. (a -> b) -> a -> b
$ Complex Algebraic
c) Unipol Rational
f Complex Algebraic -> Complex Algebraic -> Bool
forall a. Eq a => a -> a -> Bool
== Algebraic -> Algebraic -> Complex Algebraic
forall a. a -> a -> Complex a
Complex Algebraic
0 Algebraic
0
      ]

realPartPoly :: Unipol Rational -> Unipol Rational
realPartPoly :: Unipol Rational -> Unipol Rational
realPartPoly Unipol Rational
f =
  Unipol Rational -> Unipol Rational
forall r. (Eq r, Euclidean r, Division r) => Unipol r -> Unipol r
normalize (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ (Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational
2 Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Multiplicative r => r -> r -> r
* Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0) (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Unipol Rational -> Unipol Rational
rootSumPoly Unipol Rational
f Unipol Rational
f

imagPartPoly :: Unipol Rational -> Unipol Rational
imagPartPoly :: Unipol Rational -> Unipol Rational
imagPartPoly Unipol Rational
f =
  let bi :: Unipol Rational
bi =
        (Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational
2 Unipol Rational -> Unipol Rational -> Unipol Rational
forall r. Multiplicative r => r -> r -> r
* IsLabel "x" (Unipol Rational)
Unipol Rational
#x) (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
          Unipol Rational -> Unipol Rational -> Unipol Rational
rootSumPoly Unipol Rational
f ((Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Unipol Rational
forall r. Group r => r -> r
negate (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Ordinal (Arity (Unipol Rational)) -> Unipol Rational
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol Rational))
0) Unipol Rational
f)
   in (OrderedMonomial
   (MOrder (Unipol Rational)) (Arity (Unipol Rational))
 -> OrderedMonomial
      (MOrder (Unipol Rational)) (Arity (Unipol Rational)))
-> Unipol Rational -> Unipol Rational
forall poly.
IsOrderedPolynomial poly =>
(OrderedMonomial (MOrder poly) (Arity poly)
 -> OrderedMonomial (MOrder poly) (Arity poly))
-> poly -> poly
mapMonomialMonotonic ((Monomial 1 -> Identity (Monomial 1))
-> OrderedMonomial Grevlex 1
-> Identity (OrderedMonomial Grevlex 1)
forall s t. Rewrapping s t => Iso s t (Unwrapped s) (Unwrapped t)
_Wrapped ((Monomial 1 -> Identity (Monomial 1))
 -> OrderedMonomial Grevlex 1
 -> Identity (OrderedMonomial Grevlex 1))
-> ((Int -> Identity Int) -> Monomial 1 -> Identity (Monomial 1))
-> (Int -> Identity Int)
-> OrderedMonomial Grevlex 1
-> Identity (OrderedMonomial Grevlex 1)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Index (Monomial 1)
-> Traversal' (Monomial 1) (IxValue (Monomial 1))
forall m. Ixed m => Index m -> Traversal' m (IxValue m)
ix Index (Monomial 1)
0 ((Int -> Identity Int)
 -> OrderedMonomial Grevlex 1
 -> Identity (OrderedMonomial Grevlex 1))
-> Int -> OrderedMonomial Grevlex 1 -> OrderedMonomial Grevlex 1
forall a s t. Num a => ASetter s t a a -> a -> s -> t
*~ Int
2) (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
        (Ordinal (Arity (Unipol Rational)) -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall poly alg.
(IsPolynomial poly, Module (Scalar (Coefficient poly)) alg,
 Ring alg, Commutative alg) =>
(Ordinal (Arity poly) -> alg) -> poly -> alg
liftMap (Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. a -> b -> a
const (Unipol Rational -> Ordinal 1 -> Unipol Rational)
-> Unipol Rational -> Ordinal 1 -> Unipol Rational
forall a b. (a -> b) -> a -> b
$ Unipol Rational -> Unipol Rational
forall r. Group r => r -> r
negate IsLabel "x" (Unipol Rational)
Unipol Rational
#x) (Unipol Rational -> Unipol Rational)
-> Unipol Rational -> Unipol Rational
forall a b. (a -> b) -> a -> b
$
          Unipol Rational -> Unipol Rational -> Unipol Rational
rootMultPoly Unipol Rational
bi Unipol Rational
bi

-- | Choose representative element of the given interval.
representative :: (Additive r, Division r, Num r) => Interval r -> r
representative :: Interval r -> r
representative (Interval r
l r
r) = (r
l r -> r -> r
forall r. Additive r => r -> r -> r
+ r
r) r -> r -> r
forall r. Division r => r -> r -> r
/ r
2

{- | @'approximate' eps r@ returns rational number @r'@ close to @r@,
   with @abs (r - r') < eps@.
-}
approximate :: Rational -> Algebraic -> Rational
approximate :: Rational -> Algebraic -> Rational
approximate Rational
_ (Rational Rational
a) = Rational
a
approximate Rational
err (Algebraic Unipol Rational
_ [Unipol Rational]
ss Interval Rational
int) =
  Interval Rational -> Rational
forall r. (Additive r, Division r, Num r) => Interval r -> r
representative (Interval Rational -> Rational) -> Interval Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Rational
-> [Unipol Rational] -> Interval Rational -> Interval Rational
iterateImproveStrum Rational
err [Unipol Rational]
ss Interval Rational
int

-- | Same as @'approximate'@, but returns @'Fractional'@ value instead.
approxFractional :: Fractional r => Rational -> Algebraic -> r
approxFractional :: Rational -> Algebraic -> r
approxFractional Rational
eps = Rational -> r
forall a. Fractional a => Rational -> a
fromFraction (Rational -> r) -> (Algebraic -> Rational) -> Algebraic -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Algebraic -> Rational
approximate Rational
eps