module Resultant where

import Numeric.LinearAlgebra

type Poly = [Double]
  -- 多項式はこれで表す
  -- 主係数から定数までを，次数の大きい項から順に全て列挙
  -- 例 : x^2 - 2 -> [1, 0, (-2)]

deg :: Poly -> Int -- 次数計算
deg f = length f - 1


mj :: Poly -> Poly -> Int -> Matrix Double -- M_j (f, g) の計算
mj f g j = (k >< k) -- k * k 行列
           $ concat -- 行ごとの配列になっているものをくっつけて1つの配列に
           $  map (makeMsRow f k) [0..(n - j - 1)] -- 上側
           ++ map (makeMsRow g k) [0..(m - j - 1)] -- 下側
  where
    (m, n) = (deg f, deg g)
    k = m + n - 2 * j


psc :: Poly -> Poly -> Int -> Double -- PSC_j (f, g) の計算
psc f g = det . mj f g


degGcd :: Poly -> Poly -> Int
degGcd f g = head $ filter (\j -> psc f g j /= 0) [0..] -- PSC が最初に 0 じゃなくなる j


makeMsRow :: Poly -> Int -> Int -> Poly -- M_j の各行を生成
makeMsRow f k = flip add0last k . add0head f -- 右に n ずらしてから長さ k の配列を作る


add0head :: Poly -> Int -> Poly
add0head a 0 = a
add0head a n
  | n < 0     = undefined
  | otherwise = 0 : add0head a (n - 1) -- n 個の 0 を先頭に加える


add0last :: Poly -> Int -> Poly
add0last a n
  | dif == 0  = a
  | dif <  0  = add0last (init a) n -- 長い分を削る
  | otherwise = a ++ zeros dif      -- 短い分 0 を後ろに加える
  where
    dif = n - length a
    zeros :: Int -> Poly -- n 個の 0 のリストを生成
    zeros 0 = []
    zeros n
      | n < 0     = undefined
      | otherwise = 0 : zeros (n - 1)