-- |
-- Module: Symtegration.Integration.Monomial
-- Description: Integration of functions in a monomial extension.
-- Copyright: Copyright 2025 Yoo Chung
-- License: Apache-2.0
-- Maintainer: dev@chungyc.org
--
-- Supports the integration of functions in a monomial extension.
-- Specifically, for a differential field \(k\) with a monomial
-- extension \(k(t)\), attempts to derive an integral for functions
-- of the form
--
-- \[ f(t) = \frac{a(t)}{d(t)} \]
--
-- where \(t\) is a monomial over \(k\), and \(a\) and \(f\) are polynomials of \(t\).
-- Here, a monomial \(t\) with respect to a derivation \(D\) refers to a transcendental
-- over \(k\) for which \(Dt\) is a polynomial of \(t\).  For example, \(t = \tan x\)
-- is a monomial over the real numbers with respect to \(D=\frac{d}{dx}\),
-- since \(\tan x\) is transcendental and \(Dt = t^2 + 1\).
-- It does not mean a polynomial with only one term such as \(c x^n\).
--
-- A function is /simple/ with respect to a derivation \(D\) if it has
-- a denominator \(d\) which is normal, i.e., \(\gcd(d, Dd) = 1\).
-- A function is /reduced/ with respect to a derivation \(D\) if it has
-- a denominator \(d\) which is special, i.e., \(\gcd(d, Dd) = d\).
module Symtegration.Integration.Monomial
  ( hermiteReduce,
    polynomialReduce,
    residueReduce,
  )
where

import Data.List (find)
import Data.Monoid (Sum (..))
import Symtegration.Polynomial
import Symtegration.Polynomial.Differential
import Symtegration.Polynomial.Rational

-- $setup
-- >>> import Symtegration.Polynomial
-- >>> import Symtegration.Polynomial.Differential
-- >>> import Symtegration.Polynomial.Indexed
-- >>> import Symtegration.Polynomial.Rational
-- >>> import Symtegration.Symbolic.Haskell

-- | Hermite reduction for a function in a monomial extension.
--
-- Specifically, for derivation \(D\) and function \(f\),
-- returns \(g\), \(h\), and \(r\) such that
--
-- \[ f = Dg + h + r \]
--
-- where \(h\) is simple and \(r\) is reduced.
--
-- For example, with derivation \(D\) and monomial \(t\) such that \(Dt = t^2 + 1\),
-- it is the case that \(\frac{t^4 + t + 1}{t^2} = Dg + h + r\),
-- where \(g = -\frac{1}{t}\), \(h = \frac{1}{t}\), and \(r = t^2 - 1\).
--
-- >>> let deriv = extend (const 0) (power 2 + 1 :: IndexedPolynomial)
-- >>> hermiteReduce deriv $ fromPolynomials (power 4 + power 1 + 1) (power 2)
-- ([Function ((-1)) (x)],Function (1) (x),Function (x^2 + (-1)) (1))
--
-- Since it is the case that \(\frac{dt}{dx} = t^2 + 1\) for \(t = \tan x\),
-- this implies that
--
-- \[
--   \int \frac{\tan^4 x + \tan x + 1}{\tan^2 x} \, dx =
--   -\frac{1}{\tan x} + \int \left( \frac{1}{\tan x} + \tan^2 x - 1\right) \, dx
-- \]
--
-- \(g\) is returned as a list of functions which sum to \(g\)
-- instead of a single function, because the former could sometimes
-- be simpler to read.
hermiteReduce ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
  -- | Derivation \(D\).
  (p e c -> p e c) ->
  -- | Function \(f\).
  Function (p e c) ->
  -- | Hermite reduction \((g, h, r)\) of \(f\).
  ([Function (p e c)], Function (p e c), Function (p e c))
hermiteReduce :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
(p e c -> p e c)
-> Function (p e c)
-> ([Function (p e c)], Function (p e c), Function (p e c))
hermiteReduce p e c -> p e c
_ r :: Function (p e c)
r@(Function p e c
_ p e c
0) = ([], Function (p e c)
0, Function (p e c)
r)
hermiteReduce p e c -> p e c
derivation Function (p e c)
f = ([Function (p e c)]
g, Function (p e c)
h, p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
p) Function (p e c) -> Function (p e c) -> Function (p e c)
forall a. Num a => a -> a -> a
+ Function (p e c)
s)
  where
    (p e c
p, Function (p e c)
n, Function (p e c)
s) = (p e c -> p e c)
-> Function (p e c) -> (p e c, Function (p e c), Function (p e c))
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
(p e c -> p e c)
-> Function (p e c) -> (p e c, Function (p e c), Function (p e c))
canonical p e c -> p e c
derivation Function (p e c)
f
    ([Function (p e c)]
g, Function (p e c)
h, p e c
q) = (p e c -> p e c)
-> Function (p e c)
-> ([Function (p e c)], Function (p e c), p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
(p e c -> p e c)
-> Function (p e c)
-> ([Function (p e c)], Function (p e c), p e c)
hermiteReduce' p e c -> p e c
derivation Function (p e c)
n

-- | Hermite reduction on the normal part of a canonical representation.
-- This does the main work for 'hermiteReduce'.
hermiteReduce' ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
  (p e c -> p e c) ->
  Function (p e c) ->
  ([Function (p e c)], Function (p e c), p e c)
hermiteReduce' :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
(p e c -> p e c)
-> Function (p e c)
-> ([Function (p e c)], Function (p e c), p e c)
hermiteReduce' p e c -> p e c
derivation f :: Function (p e c)
f@(Function p e c
x p e c
y)
  | (Just ([Function (p e c)]
g, (p e c
a, p e c
d))) <- p e c
-> [Function (p e c)]
-> p e c
-> Maybe ([Function (p e c)], (p e c, p e c))
reduce p e c
x [] p e c
common =
      let (p e c
q, p e c
r) = p e c
a p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
d
       in ([Function (p e c)]
g, p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials p e c
r p e c
d, p e c
q)
  | Bool
otherwise = ([], Function (p e c)
f, p e c
0) -- Should not be possible.
  where
    common :: p e c
common = p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
y (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
derivation p e c
y
    (p e c
divisor, p e c
_) = p e c
y p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
common
    reduce :: p e c
-> [Function (p e c)]
-> p e c
-> Maybe ([Function (p e c)], (p e c, p e c))
reduce p e c
a [Function (p e c)]
g p e c
d
      | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
d e -> e -> Bool
forall a. Ord a => a -> a -> Bool
> e
0 = do
          let d' :: p e c
d' = p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
d (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
derivation p e c
d
          let (p e c
d'', p e c
_) = p e c
d p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
d'
          let (p e c
d''', p e c
_) = (p e c
divisor p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c -> p e c
derivation p e c
d) p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
d
          (p e c
b, p e c
c) <- p e c -> p e c -> p e c -> Maybe (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c -> Maybe (p e c, p e c)
diophantineEuclidean (-p e c
d''') p e c
d'' p e c
a
          let (p e c
b', p e c
_) = (p e c -> p e c
derivation p e c
b p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
divisor) p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
d''
          let a' :: p e c
a' = p e c
c p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c
b'
          let g' :: [Function (p e c)]
g' = p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials p e c
b p e c
d Function (p e c) -> [Function (p e c)] -> [Function (p e c)]
forall a. a -> [a] -> [a]
: [Function (p e c)]
g
          p e c
-> [Function (p e c)]
-> p e c
-> Maybe ([Function (p e c)], (p e c, p e c))
reduce p e c
a' [Function (p e c)]
g' p e c
d'
      | Bool
otherwise = ([Function (p e c)], (p e c, p e c))
-> Maybe ([Function (p e c)], (p e c, p e c))
forall a. a -> Maybe a
Just ([Function (p e c)]
g, (p e c
a, p e c
divisor))

-- | Polynomial reduction in a monomial extension with the given derivation and polynomial.
-- Specifically, for derivation \(D\) and polynomial \(p\),
-- returns polynomials \(g\) and \(h\) such that
--
-- \[ p = Dg + h \]
--
-- where the degree of \(h\) is less than the degree of \(Dt\),
-- if the latter is larger than 0.
--
-- For example, with derivation \(D\) such that \(Dt = t^2 + 1\),
-- it is the case that \(3t^4 + 3t^3 + 4 = D\left( t^3 + \frac{3}{2}t^2 - 3t \right) - 3t + 7\).
--
-- >>> let derivation = extend (const 0) (power 2 + 1 :: IndexedPolynomial)
-- >>> polynomialReduce derivation $ 3 * power 4 + 3 * power 3 + 4
-- (x^3 + (3 % 2)x^2 + (-3)x,(-3)x + 7)
polynomialReduce ::
  (Polynomial p e c, Num (p e c), Fractional c) =>
  -- | Derivation \(D\).
  (p e c -> p e c) ->
  -- | Polynomial \(p\).
  p e c ->
  -- | Polynomial reduction \((g, h)\).
  (p e c, p e c)
polynomialReduce :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
(p e c -> p e c) -> p e c -> (p e c, p e c)
polynomialReduce p e c -> p e c
derivation p e c
p
  | e
delta e -> e -> Bool
forall a. Eq a => a -> a -> Bool
== e
0 = (p e c
0, p e c
p)
  | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> Bool
forall a. Ord a => a -> a -> Bool
< e
delta = (p e c
0, p e c
p)
  | Bool
otherwise = (p e c
q0 p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
q, p e c
r)
  where
    delta :: e
delta = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree (p e c -> e) -> p e c -> e
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
derivation (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
1
    lambda :: c
lambda = p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient (p e c -> c) -> p e c -> c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
derivation (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
1
    m :: e
m = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> e
forall a. Num a => a -> a -> a
- e
delta e -> e -> e
forall a. Num a => a -> a -> a
+ e
1
    q0 :: p e c
q0 = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
p c -> c -> c
forall a. Fractional a => a -> a -> a
/ (e -> c
forall a b. (Integral a, Num b) => a -> b
fromIntegral e
m c -> c -> c
forall a. Num a => a -> a -> a
* c
lambda)) (e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
m)

    -- Explicitly delete leading terms so that they are canceled out even for symbolic coefficients.
    (p e c
q, p e c
r) = (p e c -> p e c) -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
(p e c -> p e c) -> p e c -> (p e c, p e c)
polynomialReduce p e c -> p e c
derivation (p e c -> (p e c, p e c)) -> p e c -> (p e c, p e c)
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm p e c
p p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm (p e c -> p e c
derivation p e c
q0)

-- | Resultant reduction of a simple function in a monomial extension.
--
-- Specifically, for a derivation \(D\) on \(k(t)\) and a simple function \(f \in k(t)\),
-- returns logarithmic terms \((s_i, S_i)\) such that
--
-- \[ g = \sum_i \sum_{\alpha \mid s_i(\alpha) = 0} \alpha \log S_i(\alpha, t) \]
--
-- and also whether \(f\) has an elemantary integral.  If \(f\) has an elementary integral,
-- then there is a polynomial \(h\) of the monomial \(t\) such that \(f - Dg = h\).
--
-- For example, for derivation \(D=\frac{d}{dx}\) and simple function \(f = \frac{t^2+2t+1}{t^2-t+1}\),
--
-- >>> residueReduce differentiate $ fromPolynomials (power 2 + 2 * power 1 + 1) (2 * power 2 - 2 * power 1 - 1 :: IndexedPolynomial)
-- Just ([(x^2 + ((-3) % 2)x + ((-3) % 16),[(0,(4 % 9)x + (1 % 3)),(1,((-8) % 9)x + (2 % 3))])],True)
--
-- so that
--
-- \[ g = \sum_{\alpha \mid \alpha^2-\frac{3}{2}\alpha-\frac{3}{16} = 0} \log \left( (\frac{2}{3}-\frac{8}{9}\alpha)t + \frac{4}{9}\alpha+\frac{1}{3} \right) \]
--
-- and \(\int f \, dx\) has an elementary integral, which happens to be
--
-- \[ \int f \, dx = g + \int h \, dx \]
--
-- where \(h = f-Dg\) is a polynomial of the monomial \(t\).
residueReduce ::
  ( Polynomial p e c,
    Eq (p e c),
    Num (p e c),
    Polynomial p e (p e c),
    Eq (p e (p e c)),
    Num (p e (p e c)),
    Polynomial p e (Function (p e c)),
    Eq (p e (Function (p e c))),
    Num (p e (Function (p e c))),
    Eq c,
    Fractional c
  ) =>
  -- | Derivation \(D\).
  (p e c -> p e c) ->
  -- | Function \(f\).
  Function (p e c) ->
  -- | Logarithmic terms \((s_i, S_i)\) and whether an elementary integral for \(\int f \, dx\) exists.
  Maybe ([(p e c, p e (p e c))], Bool)
residueReduce :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Polynomial p e (p e c),
 Eq (p e (p e c)), Num (p e (p e c)),
 Polynomial p e (Function (p e c)), Eq (p e (Function (p e c))),
 Num (p e (Function (p e c))), Eq c, Fractional c) =>
(p e c -> p e c)
-> Function (p e c) -> Maybe ([(p e c, p e (p e c))], Bool)
residueReduce p e c -> p e c
derivation (Function p e c
e p e c
d) = do
  -- Polynomials of monomial t.
  let (p e c
_, p e c
a) = p e c
e p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
d

  -- Coefficients are polynomials of z in rational function form.
  let d' :: p e (Function (p e c))
d' = (c -> Function (p e c)) -> p e c -> p e (Function (p e c))
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (\c
c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c)) -> p e c -> Function (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c p e c
1) p e c
d
  let a' :: p e (Function (p e c))
a' = (c -> Function (p e c)) -> p e c -> p e (Function (p e c))
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (\c
c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c)) -> p e c -> Function (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c p e c
1) p e c
a
  let zd' :: p e (Function (p e c))
zd' = (c -> Function (p e c)) -> p e c -> p e (Function (p e c))
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (\c
c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c)) -> p e c -> Function (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
1) (p e c -> p e c
derivation p e c
d)

  let (Function (p e c)
r', [p e (Function (p e c))]
rs')
        | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree (p e c -> p e c
derivation p e c
d) e -> e -> Bool
forall a. Ord a => a -> a -> Bool
<= p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
d = p e (Function (p e c))
-> p e (Function (p e c))
-> (Function (p e c), [p e (Function (p e c))])
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
p e c -> p e c -> (c, [p e c])
subresultant p e (Function (p e c))
d' (p e (Function (p e c))
a' p e (Function (p e c))
-> p e (Function (p e c)) -> p e (Function (p e c))
forall a. Num a => a -> a -> a
- p e (Function (p e c))
zd')
        | Bool
otherwise = p e (Function (p e c))
-> p e (Function (p e c))
-> (Function (p e c), [p e (Function (p e c))])
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
p e c -> p e c -> (c, [p e c])
subresultant (p e (Function (p e c))
a' p e (Function (p e c))
-> p e (Function (p e c)) -> p e (Function (p e c))
forall a. Num a => a -> a -> a
- p e (Function (p e c))
zd') p e (Function (p e c))
d'

  -- Polynomial of z.
  p e c
r <- Function (p e c) -> Maybe (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
Function (p e c) -> Maybe (p e c)
toPolynomial Function (p e c)
r'

  -- Polynomials of t with coefficients which are polynomials of z.
  [p e (p e c)]
rs <- (p e (Function (p e c)) -> Maybe (p e (p e c)))
-> [p e (Function (p e c))] -> Maybe [p e (p e c)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM ((Function (p e c) -> Maybe (p e c))
-> p e (Function (p e c)) -> Maybe (p e (p e c))
forall (p :: * -> * -> *) e c c' (m :: * -> *).
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c'),
 Monad m) =>
(c -> m c') -> p e c -> m (p e c')
mapCoefficientsM Function (p e c) -> Maybe (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
Function (p e c) -> Maybe (p e c)
toPolynomial) [p e (Function (p e c))]
rs'

  -- Splitting factorization with respect to coefficient lifting.
  -- Polynomials of z with coefficients which are polynomials of t.
  let kderiv :: p e c -> p e c
kderiv = Sum (p e c) -> p e c
forall a. Sum a -> a
getSum (Sum (p e c) -> p e c) -> (p e c -> Sum (p e c)) -> p e c -> p e c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (e -> c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall m. Monoid m => (e -> c -> m) -> p e c -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\e
ex c
c -> p e c -> Sum (p e c)
forall a. a -> Sum a
Sum (p e c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall a b. (a -> b) -> a -> b
$ p e c -> p e c
derivation (c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c p e c
1) p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
ex)
  let factors :: [(p e c, p e c)]
factors = (p e c -> p e c) -> p e c -> [(p e c, p e c)]
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
(p e c -> p e c) -> p e c -> [(p e c, p e c)]
splitSquarefreeFactor p e c -> p e c
kderiv p e c
r

  -- Whether there exists an elementary integral for the function of t in field c.
  let elementary :: Bool
elementary = ((p e c, p e c) -> Bool) -> [(p e c, p e c)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (e -> e -> Bool
forall a. Eq a => a -> a -> Bool
(==) e
0 (e -> Bool) -> ((p e c, p e c) -> e) -> (p e c, p e c) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree (p e c -> e) -> ((p e c, p e c) -> p e c) -> (p e c, p e c) -> e
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (p e c, p e c) -> p e c
forall a b. (a, b) -> a
fst) [(p e c, p e c)]
factors

  -- Derive g.  For each term (s, S), each s is a polynomial of z,
  -- and each S is a polynomial of t with coefficients which are polynomials of z.
  let specials :: [p e c]
specials = ((p e c, p e c) -> p e c) -> [(p e c, p e c)] -> [p e c]
forall a b. (a -> b) -> [a] -> [b]
map (p e c, p e c) -> p e c
forall a b. (a, b) -> b
snd [(p e c, p e c)]
factors
  [(p e c, p e (p e c))]
terms' <- ((e, p e c) -> Maybe (p e c, p e (p e c)))
-> [(e, p e c)] -> Maybe [(p e c, p e (p e c))]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM ([p e (p e c)]
-> [p e c] -> (e, p e c) -> Maybe (p e c, p e (p e c))
forall {t :: * -> *} {p :: * -> * -> *} {e} {c}.
(Foldable t, Polynomial p e c, Fractional c, Eq c) =>
t (p e (p e c))
-> [p e c] -> (e, p e c) -> Maybe (p e c, p e (p e c))
toTerm [p e (p e c)]
rs [p e c]
specials) ([(e, p e c)] -> Maybe [(p e c, p e (p e c))])
-> [(e, p e c)] -> Maybe [(p e c, p e (p e c))]
forall a b. (a -> b) -> a -> b
$ [e] -> [p e c] -> [(e, p e c)]
forall a b. [a] -> [b] -> [(a, b)]
zip [e
1 ..] [p e c]
specials

  -- Remove terms which will be equal to zero.
  let terms :: [(p e c, p e (p e c))]
terms = ((p e c, p e (p e c)) -> Bool)
-> [(p e c, p e (p e c))] -> [(p e c, p e (p e c))]
forall a. (a -> Bool) -> [a] -> [a]
filter ((p e (p e c) -> p e (p e c) -> Bool
forall a. Eq a => a -> a -> Bool
/= p e (p e c)
1) (p e (p e c) -> Bool)
-> ((p e c, p e (p e c)) -> p e (p e c))
-> (p e c, p e (p e c))
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (p e c, p e (p e c)) -> p e (p e c)
forall a b. (a, b) -> b
snd) [(p e c, p e (p e c))]
terms'

  ([(p e c, p e (p e c))], Bool)
-> Maybe ([(p e c, p e (p e c))], Bool)
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return ([(p e c, p e (p e c))]
terms, Bool
elementary)
  where
    -- Return the logarithmic term corresponding to the given special factor.
    toTerm :: t (p e (p e c))
-> [p e c] -> (e, p e c) -> Maybe (p e c, p e (p e c))
toTerm t (p e (p e c))
prs [p e c]
specials (e
i, p e c
s)
      | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
s e -> e -> Bool
forall a. Eq a => a -> a -> Bool
== e
0 = (p e c, p e (p e c)) -> Maybe (p e c, p e (p e c))
forall a. a -> Maybe a
Just (p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic p e c
s, p e (p e c)
1)
      | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
d e -> e -> Bool
forall a. Eq a => a -> a -> Bool
== e
i = (p e c, p e (p e c)) -> Maybe (p e c, p e (p e c))
forall a. a -> Maybe a
Just (p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic p e c
s, (c -> p e c) -> p e c -> p e (p e c)
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
`scale` p e c
1) p e c
d)
      | Just p e (p e c)
pr <- (p e (p e c) -> Bool) -> t (p e (p e c)) -> Maybe (p e (p e c))
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (e -> e -> Bool
forall a. Eq a => a -> a -> Bool
(==) e
i (e -> Bool) -> (p e (p e c) -> e) -> p e (p e c) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. p e (p e c) -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree) t (p e (p e c))
prs = p e c -> p e (p e c) -> [p e c] -> Maybe (p e c, p e (p e c))
forall {p :: * -> * -> *} {e} {c} {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Polynomial p e c, Polynomial p e (p e c),
 Polynomial p e (Function (p e c)), Fractional c, Fractional c,
 Eq c, Eq c, Eq (p e c), Eq (p e (Function (p e c))),
 Num (p e (p e c)), Num (p e (Function (p e c)))) =>
p e c -> p e (p e c) -> [p e c] -> Maybe (p e c, p e (p e c))
derive p e c
s p e (p e c)
pr [p e c]
specials
      | Bool
otherwise = Maybe (p e c, p e (p e c))
forall a. Maybe a
Nothing -- Should not be possible.

    -- Derive the logarithmic term given the polynomial remaindder.
    derive :: p e c -> p e (p e c) -> [p e c] -> Maybe (p e c, p e (p e c))
derive p e c
s p e (p e c)
pr [p e c]
specials = do
      -- Switch to a polynomial of z with coefficients which are polynomials of t.
      let pr' :: p e (Function (p e c))
pr' = (p e c -> Function (p e c))
-> p e (p e c) -> p e (Function (p e c))
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e (p e c) -> p e (Function (p e c)))
-> p e (p e c) -> p e (Function (p e c))
forall a b. (a -> b) -> a -> b
$ p e (p e c) -> p e (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e (p e c), Polynomial p e c, Num (p e (p e c))) =>
p e (p e c) -> p e (p e c)
switchVars p e (p e c)
pr

      let (p e (Function (p e c))
logArg', p e (Function (p e c))
_) = p e (Function (p e c))
pr' p e (Function (p e c))
-> p e (Function (p e c))
-> (p e (Function (p e c)), p e (Function (p e c)))
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e (Function (p e c))
divisor
      p e (p e c)
logArg <- (Function (p e c) -> Maybe (p e c))
-> p e (Function (p e c)) -> Maybe (p e (p e c))
forall (p :: * -> * -> *) e c c' (m :: * -> *).
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c'),
 Monad m) =>
(c -> m c') -> p e c -> m (p e c')
mapCoefficientsM Function (p e c) -> Maybe (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
Function (p e c) -> Maybe (p e c)
toPolynomial p e (Function (p e c))
logArg'

      -- Return the logarithmic term,
      -- with the argument switched back to a polynomial of t with coefficients of z.
      (p e c, p e (p e c)) -> Maybe (p e c, p e (p e c))
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic p e c
s, p e (p e c) -> p e (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e (p e c), Polynomial p e c, Num (p e (p e c))) =>
p e (p e c) -> p e (p e c)
switchVars p e (p e c)
logArg)
      where
        -- Polynomials of z.
        factors :: [p e c]
factors = p e c -> [p e c]
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
p e c -> [p e c]
squarefree (p e c -> [p e c]) -> p e c -> [p e c]
forall a b. (a -> b) -> a -> b
$ p e (p e c) -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e (p e c)
pr

        -- Polynomial of z with coefficients which are polynomials of t.
        divisor :: p e (Function (p e c))
divisor = (c -> Function (p e c)) -> p e c -> p e (Function (p e c))
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c))
-> (c -> p e c) -> c -> Function (p e c)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (c -> p e c -> p e c) -> p e c -> c -> p e c
forall a b c. (a -> b -> c) -> b -> a -> c
flip c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale p e c
1) p e c
divisor'
          where
            divisor' :: p e c
divisor' = [p e c] -> p e c
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([p e c] -> p e c) -> [p e c] -> p e c
forall a b. (a -> b) -> a -> b
$ ((Int, p e c, p e c) -> p e c) -> [(Int, p e c, p e c)] -> [p e c]
forall a b. (a -> b) -> [a] -> [b]
map (Int, p e c, p e c) -> p e c
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Fractional c, Eq (p e c), Num (p e c)) =>
(Int, p e c, p e c) -> p e c
toDivisor ([(Int, p e c, p e c)] -> [p e c])
-> [(Int, p e c, p e c)] -> [p e c]
forall a b. (a -> b) -> a -> b
$ [Int] -> [p e c] -> [p e c] -> [(Int, p e c, p e c)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [Int
1 ..] [p e c]
factors [p e c]
specials
            toDivisor :: (Int, p e c, p e c) -> p e c
toDivisor (Int
j, p e c
factor, p e c
special) = p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
factor p e c
special p e c -> Int -> p e c
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
j :: Int)

-- | Given a polynomial with coefficients which are polynomials of a different variable,
-- return an equivalent polynomial of the different variable with coefficients
-- which are polynomials of the original variable.
--
-- For example, a polynomial such as
--
-- \[ (y^2+1)x^3 + (2y^2-y)x + 2y - 1 \]
--
-- would be turned into
--
-- \[ (x^3 + 2x) y^2 - (x-2)y + x^3 - 1 \]
switchVars ::
  (Polynomial p e (p e c), Polynomial p e c, Num (p e (p e c))) =>
  -- | Polynomial of \(v_1\) with coefficients which are polynomials of \(v_2\).
  p e (p e c) ->
  -- | Polynomial of \(v_2\) with coefficients which are polynomials of \(v_1\).
  p e (p e c)
switchVars :: forall (p :: * -> * -> *) e c.
(Polynomial p e (p e c), Polynomial p e c, Num (p e (p e c))) =>
p e (p e c) -> p e (p e c)
switchVars p e (p e c)
p = Sum (p e (p e c)) -> p e (p e c)
forall a. Sum a -> a
getSum (Sum (p e (p e c)) -> p e (p e c))
-> Sum (p e (p e c)) -> p e (p e c)
forall a b. (a -> b) -> a -> b
$ (e -> p e c -> Sum (p e (p e c)))
-> p e (p e c) -> Sum (p e (p e c))
forall m. Monoid m => (e -> p e c -> m) -> p e (p e c) -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\e
e p e c
c -> p e (p e c) -> Sum (p e (p e c))
forall a. a -> Sum a
Sum (p e (p e c) -> Sum (p e (p e c)))
-> p e (p e c) -> Sum (p e (p e c))
forall a b. (a -> b) -> a -> b
$ (c -> p e c) -> p e c -> p e (p e c)
forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients (\c
c' -> c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c' (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
e) p e c
c) p e (p e c)
p