module Numeric.AlgebraicReal.Class where
import Data.Ratio
import qualified Data.Vector as V

-- | 整域
--
-- \(1 \ne 0\) となる可換環で、非自明な零因子を持たないものを整域という。
class (Num a) => IntegralDomain a where
  -- | 除算
  --
  -- @a@ が @b@ の倍数であるとき、 @divide a b@ は @a = b * c@ となる @c@ を返す。
  -- @a@ が @b@ の倍数でない時の挙動は規定しない。
  divide :: a -> a -> a

infixl 7 `divide`

instance (Integral a) => IntegralDomain (Ratio a) where
  divide = (/)

instance IntegralDomain Integer where
  divide = div

{-
With UndecidableInstances:
instance (Fractional a) => IntegralDomain a where
  divide = (/)
-}

class (IntegralDomain a) => GCDDomain a where
  gcdD :: a -> a -> a
  contentV :: V.Vector a -> a
  contentV = V.foldr gcdD 0

-- | 体の場合の gcdD 関数のデフォルト実装
fieldGcd :: (Eq a, Fractional a) => a -> a -> a
fieldGcd 0 0 = 0
fieldGcd _ _ = 1

-- | 体の場合の contentV 関数のデフォルト実装
fieldContentV :: (Eq a, Fractional a) => V.Vector a -> a
fieldContentV xs | V.null xs = 0
                 | otherwise = V.last xs

instance (Integral a) => GCDDomain (Ratio a) where
  gcdD = fieldGcd
  contentV = fieldContentV

instance GCDDomain Integer where
  gcdD = gcd
  contentV xs | V.null xs = 0
              | V.last xs < 0 = negate (gcdV 0 xs) -- 内容の符号に、最高次の係数の符号を反映させる
              | otherwise = gcdV 0 xs -- 短絡評価を考えなければ foldr gcd 0 でも良い
    where
      -- foldl/foldr と gcd の組み合わせでは GCD が 1 になっても残りの部分が評価される。
      -- 列の途中で GCD が 1 になれば全体の GCD は 1 で確定なので、そういう短絡評価する。
      gcdV 1 _ = 1
      gcdV a v | V.null v = a
               | otherwise = gcdV (gcdD (V.last v) a) (V.init v)

{-
With UndecidableInstances:
instance (Eq a, Fractional a) => GCDDomain a where
  gcdD = fieldGcd
  contentV = fieldContentV
-}

-- ユークリッド整域
class (Num a) => EuclideanDomain a where
  -- a, b に対して、 (q,r) = divModD a b は a = q * b + r を満たす
  divModD :: a -> a -> (a, a)
  modD :: a -> a -> a
  modD x y = snd (divModD x y)

infixl 7 `modD`

instance EuclideanDomain Integer where
  divModD = divMod

{- あってもなくても良い：
instance (Integral a) => EuclideanDomain (Ratio a) where
  divModD x y = (x / y, 0)

With UndecidableInstances:
instance (Fractional a) => EuclideanDomain a where
  divModD x y = (x / y, 0)
-}

...
module Numeric.AlgebraicReal.UniPoly where
import Numeric.AlgebraicReal.Class
...

unscaleP :: (Eq a, IntegralDomain a) => a -> UniPoly a -> UniPoly a
unscaleP a f = mapCoeff (`divide` a) f

-- | 多項式の内容を計算する
content :: (GCDDomain a) => UniPoly a -> a
content (UniPoly xs) = contentV xs

-- | 多項式の内容と原始部分を計算する
contentAndPrimitivePart :: (Eq a, GCDDomain a) => UniPoly a -> (a, UniPoly a)
contentAndPrimitivePart f
  | c == 1 = (c, f)
  | otherwise = (c, unscaleP c f)
  where c = content f

-- | 多項式の原始部分を計算する
primitivePart :: (Eq a, GCDDomain a) => UniPoly a -> UniPoly a
primitivePart = snd . contentAndPrimitivePart

instance (Eq a, IntegralDomain a) => IntegralDomain (UniPoly a) where
  divide f g
    | g == 0 = error "divide: divide by zero"
    | degree f < degree g = zeroP -- f should be zero
    | otherwise = loop zeroP f
    where
      l = degree' f - degree' g + 1
      b = leadingCoefficient g
      -- invariant: f == q * g + r
      loop q r | degree r < degree g = q -- r should be zero
               | otherwise = loop (q + q') (r - q' * g)
        where q' = unscaleP b (UniPoly (V.drop (degree' g) (coeff r)))

gcd_subresultantPRS :: (Eq a, IntegralDomain a) => UniPoly a -> UniPoly a -> UniPoly a
gcd_subresultantPRS f 0 = f
gcd_subresultantPRS f g
  | degree f < degree g = gcd_subresultantPRS g f
  | otherwise = case pseudoModP f g of
      0 -> g
      rem -> let !d = degree' f - degree' g
                 !s = (-1)^(d + 1) * rem
             in loop d (-1) g s
  where
    loop !_ _ f 0 = f
    loop !d psi f g = case pseudoModP f g of
      0 -> g
      rem -> let !d' = degree' f - degree' g
                 !c = leadingCoefficient f
                 !psi' | d == 0 = psi
                       | d > 0 = ((-c)^d) `divide` (psi^(d-1))
                 !beta = -c * psi' ^ d'
                 !s = unscaleP beta rem
             in loop d' psi' g s

instance (Eq a, GCDDomain a) => GCDDomain (UniPoly a) where
  gcdD x y = let (xc,xp) = contentAndPrimitivePart x
                 (yc,yp) = contentAndPrimitivePart y
             in scaleP (gcdD xc yc) $ primitivePart (gcd_subresultantPRS xp yp)
  contentV xs | V.null xs = 0
              | otherwise = scaleP (contentV (V.map fst ys)) $ gcdV 0 (V.map snd ys)
    where
      ys = V.map contentAndPrimitivePart xs
      gcdV 1 _ = 1
      gcdV a v | V.null v = a
               | otherwise = gcdV (primitivePart $ gcd_subresultantPRS (V.last v) a) (V.init v)

instance (Eq a, Fractional a, GCDDomain a) => EuclideanDomain (UniPoly a) where
  divModD = divModP

-- Resultant.hs に追加

-- subresultantPRS' f g = (b,r) : subresultantPRS' g r
-- where lc(g)^l * f = q * g + b * r, l = deg f - deg g + 1
subresultantPRS' :: (Eq a, IntegralDomain a) => UniPoly a -> UniPoly a -> [(a,UniPoly a)]
subresultantPRS' _ 0 = []
subresultantPRS' f g = case pseudoModP f g of
  0 -> []
  rem -> let !d = degree' f - degree' g
             !s = (-1)^(d + 1) * rem
         in ((-1)^(d + 1), s) : loop d (-1) g s
  where
    loop !_ _ _ 0 = []
    loop d psi f g = case pseudoModP f g of
      0 -> []
      rem -> let !d' = degree' f - degree' g
                 !c = leadingCoefficient f
                 !psi' = (-c)^d `divide` psi^(d-1)
                 !beta = -c * psi' ^ d'
                 !s = unscaleP beta rem
             in (beta,s) : loop d' psi' g s

-- AlgReal.hs の negativePRS を置き換え

-- | Negative polynomial remainder sequence using subresultant PRS
--
-- Assumption: 'degree f > degree g'
negativePRS :: (Ord a, IntegralDomain a) => UniPoly a -> UniPoly a -> [UniPoly a]
negativePRS f g = f : g : loop 1 f 1 g (subresultantPRS' f g)
  where
    loop !_ _ !_ _ [] = []
    loop !s f !t g ((b,x):xs)
      -- b * (lc g)^(degree f - degree g + 1) * s > 0
      | sign b * (sign $ leadingCoefficient g)^(degree' f - degree' g + 1) * s > 0 = -x : loop t g (-1) x xs
      | otherwise = x : loop t g 1 x xs

squareFree :: (Eq a, GCDDomain a) => UniPoly a -> UniPoly a
squareFree f = f `divide` gcdD f (diffP f)

data AlgReal = FromRat !Rational
             | AlgReal !(UniPoly Rational) !Int !Rational !Rational
  deriving (Show)

data AlgReal = FromRat !Rational
             | AlgReal !(UniPoly Integer) !Int !Rational !Rational
  deriving (Show)

definingPolynomial :: AlgReal -> UniPoly Rational
definingPolynomial (FromRat x) = ind - constP x
definingPolynomial (AlgReal p _ _ _) = p

definingPolynomial :: AlgReal -> UniPoly Integer
definingPolynomial (FromRat x) = fromCoeff (V.fromList [- numerator x, denominator x])
definingPolynomial (AlgReal p _ _ _) = p

-- UniPoly.hs に追加

valueAtZQ :: Rational -> UniPoly Integer -> Rational
valueAtZQ t (UniPoly xs) = V.foldr' (\a b -> fromInteger a + t * b) 0 xs

-- UniPoly.hs に追加
-- homogeneousValueAt x y (a_n X^n + ... + a_1 X + a_0)
-- = (a_n x^n + a_{n-1} x^{n-1} y + ... + a_1 x y^{n-1} + a_0 y^n, y^n)
homogeneousValueAt :: (Eq a, Num a) => a -> a -> UniPoly a -> (a, a)
homogeneousValueAt num den f@(UniPoly coeff)
  | f == 0 = (0, 1)
  | otherwise = (V.foldr' (\x y -> x + num * y) 0 $ V.zipWith (*) coeff denseq, V.head denseq)
  where
    -- numseq = V.iterateN (V.length coeff) (* num) 1
    denseq = V.reverse (V.iterateN (V.length coeff) (* den) 1)

-- homogeneousCompP (a_n X^n + ... + a_1 X + a_0) g d
-- = (a_n g^n + a_{n-1} g^{n-1} d + ... + a_1 g d^{n-1} + a_0 d^n, d^n)
homogeneousCompP :: (Eq a, Num a) => UniPoly a -> UniPoly a -> a -> (UniPoly a, a)
homogeneousCompP f@(UniPoly coeff) g den
  | f == 0 = (0, 1)
  | otherwise = (V.foldr' (\x y -> constP x + g * y) 0 $ V.zipWith (*) coeff denseq, V.head denseq)
  where
    denseq = V.reverse (V.iterateN (V.length coeff) (* den) 1)

-- AlgReal.hs に追加

signAtZQ :: Rational -> UniPoly Integer -> Int
signAtZQ x p = sign (valueAtZQ x p)

isZeroAtZQ :: Rational -> UniPoly Integer -> Bool
isZeroAtZQ x p = valueAtZQ x p == 0

intervalsWithSign :: UniPoly Integer -> Int -> Rational -> Rational -> [Interval]
mkAlgReal :: UniPoly Integer -> Rational -> Rational -> AlgReal
mkAlgRealWithInterval :: UniPoly Integer -> Interval -> AlgReal
realRoots :: UniPoly Integer -> [AlgReal]
realRootsBetween :: UniPoly Integer -> ExtReal Rational -> ExtReal Rational -> [AlgReal]
mkAlgRealWithCReal :: UniPoly Integer -> CReal -> AlgReal

realRootsBetweenQ :: UniPoly Rational -> ExtReal Rational -> ExtReal Rational -> [AlgReal]
realRootsBetweenQ f lb ub
  | f == 0 = error "realRoots: zero"
  | degree' f == 0 = []
  | otherwise = let commonDenominator = V.foldl' (\a b -> lcm a (denominator b)) 1 (coeff f)
                    fz = primitivePart $ mapCoeff (\x -> numerator x * (commonDenominator `div` denominator x)) f
                in realRootsBetween fz lb ub

rootBound :: UniPoly Integer -> Rational
rootBound f | f == 0 = error "rootBound: polynomial is zero"
            | otherwise = 1 + (V.maximum $ V.map (abs . (% lc)) $ V.init $ coeff f)
  where lc = leadingCoefficient f

instance Num AlgReal where
...
  FromRat k + AlgReal f s a b = mkAlgReal (fst $ homogeneousCompP f (definingPolynomial (FromRat k)) (denominator k)) (a + k) (b + k)
...
  FromRat k - AlgReal f s a b = mkAlgReal (fst $ homogeneousCompP f (- definingPolynomial (FromRat k)) (denominator k)) (k - b) (k - a)
  AlgReal f s a b - FromRat k = mkAlgReal (fst $ homogeneousCompP f (definingPolynomial (FromRat (-k))) (denominator k)) (a - k) (b - k)
...
  FromRat k * AlgReal f s a b
    | k == 0 = 0
    | k > 0 = AlgReal f_x_k s (a * k) (b * k)
    | k < 0 = AlgReal f_x_k (-s) (b * k) (a * k)
    where f_x_k = fst $ homogeneousValueAt (scaleP (denominator k) ind) (fromInteger $ numerator k) (mapCoeff fromInteger f) -- f(x/k)
...

instance IntegralDomain AlgReal where
  divide = (/)

instance GCDDomain AlgReal where
  gcdD = fieldGcd
  contentV = fieldContentV