-- |
-- Module: Symtegration.Integration.Rational
-- Description: Integration of rational functions.
-- Copyright: Copyright 2025 Yoo Chung
-- License: Apache-2.0
-- Maintainer: dev@chungyc.org
--
-- Integrates rational functions.
-- Rational functions are ratios of two polynomials, not functions of rational numbers.
-- Only rational number coefficients are supported.
module Symtegration.Integration.Rational
  ( -- * Integration
    integrate,

    -- * Algorithms

    -- | Algorithms used for integrating rational functions.
    hermiteReduce,
    rationalIntegralLogTerms,
    complexLogTermToAtan,
    complexLogTermToRealTerm,

    -- * Types
    RationalFunction,
  )
where

import Data.Foldable (asum)
import Data.List (find, intersect)
import Data.Monoid (Sum (..))
import Data.Text (Text)
import Symtegration.Polynomial hiding (integrate)
import Symtegration.Polynomial qualified as Polynomial
import Symtegration.Polynomial.Indexed
import Symtegration.Polynomial.Rational as Rational
import Symtegration.Polynomial.Solve
import Symtegration.Polynomial.Symbolic
import Symtegration.Symbolic
import Symtegration.Symbolic.Simplify

-- $setup
-- >>> :set -w
-- >>> import Symtegration.Polynomial hiding (integrate)
-- >>> import Symtegration.Polynomial.Indexed
-- >>> import Symtegration.Polynomial.Rational
-- >>> import Symtegration.Symbolic.Haskell
-- >>> import Symtegration.Symbolic.Simplify

-- | Integrate a ratio of two polynomials with rational number coefficients.
--
-- For example,
--
-- >>> let p = "x" ** 7 - 24 * "x" ** 4 - 4 * "x" ** 2 + 8 * "x" - 8
-- >>> let q = "x" ** 8 + 6 * "x" ** 6 + 12 * "x" ** 4 + 8 * "x" ** 2
-- >>> toHaskell . simplify <$> integrate "x" (p / q)
-- Just "3 / (2 + x ** 2) + (4 + 8 * x ** 2) / (4 * x + 4 * x ** 3 + x ** 5) + log x"
--
-- so that
--
-- \[\int \frac{x^7-24x^4-4x^2+8x-8}{x^8+6x^6+12x^4+8x^2} \, dx = \frac{3}{x^2+2} + \frac{8x^2+4}{x^5+4x^3+4x} + \log x\]
--
-- For another example,
--
-- >>> let f = 36 / ("x" ** 5 - 2 * "x" ** 4 - 2 * "x" ** 3 + 4 * "x" ** 2 + "x" - 2)
-- >>> toHaskell . simplify <$> integrate "x" f
-- Just "(-4) * log (8 + 8 * x) + 4 * log (16 + (-8) * x) + (6 + 12 * x) / ((-1) + x ** 2)"
--
-- so that
--
-- \[\int \frac{36}{x^5-2x^4-2x^3+4x^2+x-2} \, dx = \frac{12x+6}{x^2-1} + 4 \log \left( x - 2 \right) - 4 \log \left( x + 1 \right)\]
--
-- This function will attempt to find a real function integral if it can,
-- but if it cannot, it will try to find an integral which includes complex logarithms.
integrate :: Text -> Expression -> Maybe Expression
integrate :: Text -> Expression -> Maybe Expression
integrate Text
v Expression
e
  | (Expression
x :/: Expression
y) <- Expression
e',
    (Just P Int Rational
n) <- (Text -> Maybe (P Int Rational), Expression -> Maybe Rational)
-> Expression -> Maybe (P Int Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
(Text -> Maybe (p e c), Expression -> Maybe c)
-> Expression -> Maybe (p e c)
fromExpression (Text
-> (Text -> Maybe (P Int Rational), Expression -> Maybe Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
Text -> (Text -> Maybe (p e c), Expression -> Maybe c)
forVariable Text
v) Expression
x,
    (Just P Int Rational
d) <- (Text -> Maybe (P Int Rational), Expression -> Maybe Rational)
-> Expression -> Maybe (P Int Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
(Text -> Maybe (p e c), Expression -> Maybe c)
-> Expression -> Maybe (p e c)
fromExpression (Text
-> (Text -> Maybe (P Int Rational), Expression -> Maybe Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
Text -> (Text -> Maybe (p e c), Expression -> Maybe c)
forVariable Text
v) Expression
y,
    P Int Rational
d P Int Rational -> P Int Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= P Int Rational
0 =
      P Int Rational -> P Int Rational -> Maybe Expression
integrate' P Int Rational
n P Int Rational
d
  | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing
  where
    e' :: Expression
e' = Text -> Expression -> Expression
simplifyForVariable Text
v Expression
e
    integrate' :: P Int Rational -> P Int Rational -> Maybe Expression
integrate' P Int Rational
n P Int Rational
d = Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
(+) Expression
reduced (Expression -> Expression)
-> (Expression -> Expression) -> Expression -> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
(+) Expression
poly (Expression -> Expression) -> Maybe Expression -> Maybe Expression
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Maybe Expression
logs
      where
        -- Integrals directly from Hermite reduction.
        ([Function (P Int Rational)]
g, Function (P Int Rational)
h) = Function (P Int Rational)
-> ([Function (P Int Rational)], Function (P Int Rational))
hermiteReduce (Function (P Int Rational)
 -> ([Function (P Int Rational)], Function (P Int Rational)))
-> Function (P Int Rational)
-> ([Function (P Int Rational)], Function (P Int Rational))
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational -> Function (P Int Rational)
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 Int Rational
n P Int Rational
d
        reduced :: Expression
reduced = [Expression] -> Expression
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Expression] -> Expression) -> [Expression] -> Expression
forall a b. (a -> b) -> a -> b
$ (Function (P Int Rational) -> Expression)
-> [Function (P Int Rational)] -> [Expression]
forall a b. (a -> b) -> [a] -> [b]
map Function (P Int Rational) -> Expression
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Real c) =>
Function (p e c) -> Expression
fromRationalFunction [Function (P Int Rational)]
g

        -- Integrate polynomials left over from the Hermite reduction.
        Rational.Function P Int Rational
numer P Int Rational
denom = Function (P Int Rational)
h
        (P Int Rational
q, P Int Rational
r) = P Int Rational
numer P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
denom
        poly :: Expression
poly = Text -> (Rational -> Expression) -> P Int Rational -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Rational -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient (P Int Rational -> Expression) -> P Int Rational -> Expression
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
p e c -> p e c
Polynomial.integrate P Int Rational
q

        -- Derive the log terms in the integral.
        h' :: Function (P Int Rational)
h' = P Int Rational -> P Int Rational -> Function (P Int Rational)
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 Int Rational
r P Int Rational
denom
        logTerms :: Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
logTerms = Function (P Int Rational)
-> Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
rationalIntegralLogTerms Function (P Int Rational)
h'
        logs :: Maybe Expression
logs = [Maybe Expression] -> Maybe Expression
forall (t :: * -> *) (f :: * -> *) a.
(Foldable t, Alternative f) =>
t (f a) -> f a
asum [Maybe Expression
realLogs, Maybe Expression
complexLogs] :: Maybe Expression

        -- Try to integrate into real functions first.
        realLogs :: Maybe Expression
realLogs
          | (Just [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms) <- Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
logTerms = [Expression] -> Expression
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Expression] -> Expression)
-> Maybe [Expression] -> Maybe Expression
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((P Int Rational, IndexedPolynomialWith (P Int Rational))
 -> Maybe Expression)
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
-> Maybe [Expression]
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 (Text
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> Maybe Expression
complexLogTermToRealExpression Text
v) [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms
          | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing

        -- If it cannot be integrated into real functions, allow complex logarithms.
        complexLogs :: Maybe Expression
complexLogs
          | (Just [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms) <- Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
logTerms = [Expression] -> Expression
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Expression] -> Expression)
-> Maybe [Expression] -> Maybe Expression
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((P Int Rational, IndexedPolynomialWith (P Int Rational))
 -> Maybe Expression)
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
-> Maybe [Expression]
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 (Text
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> Maybe Expression
complexLogTermToComplexExpression Text
v) [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms
          | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing

        fromRationalFunction :: Function (p e c) -> Expression
fromRationalFunction (Rational.Function p e c
u p e c
w) = Expression
u' Expression -> Expression -> Expression
forall a. Fractional a => a -> a -> a
/ Expression
w'
          where
            u' :: Expression
u' = Text -> (c -> Expression) -> p e c -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v c -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient p e c
u
            w' :: Expression
w' = Text -> (c -> Expression) -> p e c -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v c -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient p e c
w

-- | Rational functions which can be integrated by this module.
type RationalFunction = Rational.Function IndexedPolynomial

-- | Applies Hermite reduction to a rational function.
-- Returns a list of rational functions whose sums add up to the integral
-- and a rational function which remains to be integrated.
-- Only rational functions with rational number coefficients and
-- where the numerator and denominator are coprime are supported.
--
-- Specifically, for rational function \(f = \frac{A}{D}\),
-- where \(A\) and \(D\) are coprime polynomials, then for return value @(gs, h)@,
-- the sum of @gs@ is equal to \(g\) and @h@ is equal to \(h\) in the following:
--
-- \[ \frac{A}{D} = \frac{dg}{dx} + h \]
--
-- This is equivalent to the following:
--
-- \[ \int \frac{A}{D} \, dx = g + \int h \, dx \]
--
-- If preconditions are satisfied, i.e., \(D \neq 0\) and \(A\) and \(D\) are coprime,
-- then \(h\) will have a squarefree denominator.
--
-- For example,
--
-- >>> let p = power 7 - 24 * power 4 - 4 * power 2 + 8 * power 1 - 8 :: IndexedPolynomial
-- >>> let q = power 8 + 6 * power 6 + 12 * power 4 + 8 * power 2 :: IndexedPolynomial
-- >>> hermiteReduce $ fromPolynomials p q
-- ([Function (3) (x^2 + 2),Function (8x^2 + 4) (x^5 + 4x^3 + 4x)],Function (1) (x))
--
-- so that
--
-- \[\int \frac{x^7-24x^4-4x^2+8x-8}{x^8+6x^6+12x^4+8x^2} \, dx = \frac{3}{x^2+2}+\frac{8x^2+4}{x^5+4x^3+4x}+\int \frac{1}{x} \, dx\]
--
-- \(g\) is returned as a list of rational functions which sum to \(g\)
-- instead of a single rational function, because the former could sometimes
-- be simpler to read.
hermiteReduce :: RationalFunction -> ([RationalFunction], RationalFunction)
hermiteReduce :: Function (P Int Rational)
-> ([Function (P Int Rational)], Function (P Int Rational))
hermiteReduce h :: Function (P Int Rational)
h@(Rational.Function P Int Rational
_ P Int Rational
0) = ([], Function (P Int Rational)
h)
hermiteReduce h :: Function (P Int Rational)
h@(Rational.Function P Int Rational
x P Int Rational
y)
  | (Just ([Function (P Int Rational)], Function (P Int Rational))
z) <- P Int Rational
-> [Function (P Int Rational)]
-> P Int Rational
-> Maybe ([Function (P Int Rational)], Function (P Int Rational))
reduce P Int Rational
x [] P Int Rational
common = ([Function (P Int Rational)], Function (P Int Rational))
z
  | Bool
otherwise = ([], Function (P Int Rational)
h) -- Should never happen, but a fallback if it does.
  where
    common :: P Int Rational
common = P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (P Int Rational -> P Int Rational)
-> P Int Rational -> P Int Rational
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational -> P Int Rational
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 Int Rational
y (P Int Rational -> P Int Rational)
-> P Int Rational -> P Int Rational
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int Rational
y
    (P Int Rational
divisor, P Int Rational
_) = P Int Rational
y P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
common
    reduce :: P Int Rational
-> [Function (P Int Rational)]
-> P Int Rational
-> Maybe ([Function (P Int Rational)], Function (P Int Rational))
reduce P Int Rational
a [Function (P Int Rational)]
g P Int Rational
d
      | P Int Rational -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree P Int Rational
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = do
          let d' :: P Int Rational
d' = P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (P Int Rational -> P Int Rational)
-> P Int Rational -> P Int Rational
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational -> P Int Rational
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 Int Rational
d (P Int Rational -> P Int Rational)
-> P Int Rational -> P Int Rational
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int Rational
d
          let (P Int Rational
d'', P Int Rational
_) = P Int Rational
d P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
d'
          let (P Int Rational
d''', P Int Rational
_) = (P Int Rational
divisor P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
* P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int Rational
d) P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
d
          (P Int Rational
b, P Int Rational
c) <- P Int Rational
-> P Int Rational
-> P Int Rational
-> Maybe (P Int Rational, P Int Rational)
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 Int Rational
d''') P Int Rational
d'' P Int Rational
a
          let (P Int Rational
b', P Int Rational
_) = (P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int Rational
b P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
* P Int Rational
divisor) P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
d''
          let a' :: P Int Rational
a' = P Int Rational
c P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
- P Int Rational
b'
          let g' :: [Function (P Int Rational)]
g' = P Int Rational -> P Int Rational -> Function (P Int Rational)
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 Int Rational
b P Int Rational
d Function (P Int Rational)
-> [Function (P Int Rational)] -> [Function (P Int Rational)]
forall a. a -> [a] -> [a]
: [Function (P Int Rational)]
g
          P Int Rational
-> [Function (P Int Rational)]
-> P Int Rational
-> Maybe ([Function (P Int Rational)], Function (P Int Rational))
reduce P Int Rational
a' [Function (P Int Rational)]
g' P Int Rational
d'
      | Bool
otherwise = ([Function (P Int Rational)], Function (P Int Rational))
-> Maybe ([Function (P Int Rational)], Function (P Int Rational))
forall a. a -> Maybe a
Just ([Function (P Int Rational)]
g, P Int Rational -> P Int Rational -> Function (P Int Rational)
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 Int Rational
a P Int Rational
divisor)

-- | For rational function \(\frac{A}{D}\), where \(\deg(A) < \deg(D)\),
-- and \(D\) is non-zero, squarefree, and coprime with \(A\),
-- returns the components which form the logarithmic terms of \(\int \frac{A}{D} \, dx\).
-- Specifically, when a list of \((Q_i(t), S_i(t, x))\) is returned,
-- where \(Q_i(t)\) are polynomials of \(t\) and \(S_i(t, x)\) are polynomials of \(x\)
-- with coefficients formed from polynomials of \(t\), then
--
-- \[
-- \int \frac{A}{D} \, dx = \sum_{i=1}^n \sum_{a \in \{t \mid Q_i(t) = 0\}} a \log \left(S_i(a,x)\right)
-- \]
--
-- For example,
--
-- >>> let p = power 4 - 3 * power 2 + 6 :: IndexedPolynomial
-- >>> let q = power 6 - 5 * power 4 + 5 * power 2 + 4 :: IndexedPolynomial
-- >>> let f = fromPolynomials p q
-- >>> let gs = rationalIntegralLogTerms f
-- >>> length <$> gs
-- Just 1
-- >>> fst . head <$> gs
-- Just x^2 + (1 % 4)
-- >>> foldTerms (\e c -> show (e, c) <> " ") . snd . head <$> gs
-- Just "(0,792x^2 + (-16)) (1,(-2440)x^3 + 32x) (2,(-400)x^2 + 7) (3,800x^3 + (-14)x) "
--
-- so it is the case that
--
-- \[
-- \int \frac{x^4-3x^2+6}{x^6-5x^4+5x^2+4} \, dx
-- = \sum_{a \mid a^2+\frac{1}{4} = 0} a \log \left( (800a^3-14a)x^3+(-400a^2+7)x^2+(-2440a^3+32a)x + 792a^2-16 \right)
-- \]
--
-- It may return 'Nothing' if \(\frac{A}{D}\) is not in the expected form.
rationalIntegralLogTerms ::
  RationalFunction ->
  Maybe [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
rationalIntegralLogTerms :: Function (P Int Rational)
-> Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
rationalIntegralLogTerms (Rational.Function P Int Rational
a P Int Rational
d) = do
  -- For A/D, get the resultant and subresultant polynomial remainder sequence
  -- for D and (A - t * D').
  let sa :: P Int (Function (P Int Rational))
sa = (Rational -> Function (P Int Rational))
-> P Int Rational -> P Int (Function (P Int Rational))
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 Rational -> Function (P Int Rational)
forall a. Fractional a => Rational -> a
fromRational P Int Rational
a
  let sd :: P Int (Function (P Int Rational))
sd = (Rational -> Function (P Int Rational))
-> P Int Rational -> P Int (Function (P Int Rational))
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 Rational -> Function (P Int Rational)
forall a. Fractional a => Rational -> a
fromRational P Int Rational
d
  let t :: Function (P Int Rational)
t = P Int Rational -> Function (P Int Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (P Int Rational -> Function (P Int Rational))
-> P Int Rational -> Function (P Int Rational)
forall a b. (a -> b) -> a -> b
$ Int -> P Int Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
  let (Function (P Int Rational)
resultant, [P Int (Function (P Int Rational))]
prs) = P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
-> (Function (P Int Rational), [P Int (Function (P Int Rational))])
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 Int (Function (P Int Rational))
sd (P Int (Function (P Int Rational))
 -> (Function (P Int Rational),
     [P Int (Function (P Int Rational))]))
-> P Int (Function (P Int Rational))
-> (Function (P Int Rational), [P Int (Function (P Int Rational))])
forall a b. (a -> b) -> a -> b
$ P Int (Function (P Int Rational))
sa P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
forall a. Num a => a -> a -> a
- Function (P Int Rational)
-> P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Function (P Int Rational)
t (P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int (Function (P Int Rational))
sd)

  -- Turn rational functions into polynomials if possible.
  -- When the preconditions are satisfied, these should all be polynomials.
  IndexedPolynomialWith (P Int Rational)
sd' <- (Function (P Int Rational) -> Maybe (P Int Rational))
-> P Int (Function (P Int Rational))
-> Maybe (IndexedPolynomialWith (P Int Rational))
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 Int Rational) -> Maybe (P Int Rational)
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 Int (Function (P Int Rational))
sd
  P Int Rational
resultant' <- Function (P Int Rational) -> Maybe (P Int Rational)
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 Int Rational)
resultant
  [IndexedPolynomialWith (P Int Rational)]
prs' <- (P Int (Function (P Int Rational))
 -> Maybe (IndexedPolynomialWith (P Int Rational)))
-> [P Int (Function (P Int Rational))]
-> Maybe [IndexedPolynomialWith (P Int Rational)]
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 Int Rational) -> Maybe (P Int Rational))
-> P Int (Function (P Int Rational))
-> Maybe (IndexedPolynomialWith (P Int Rational))
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 Int Rational) -> Maybe (P Int Rational)
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 Int (Function (P Int Rational))]
prs :: Maybe [IndexedPolynomialWith IndexedPolynomial]

  -- Derive what make up the log terms in the integral.
  let qs :: [P Int Rational]
qs = P Int Rational -> [P Int Rational]
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 Int Rational
resultant' :: [IndexedPolynomial]
  let terms :: [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms = (Int
 -> P Int Rational
 -> (P Int Rational, IndexedPolynomialWith (P Int Rational)))
-> [Int]
-> [P Int Rational]
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (IndexedPolynomialWith (P Int Rational)
-> [IndexedPolynomialWith (P Int Rational)]
-> Int
-> P Int Rational
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
toTerm IndexedPolynomialWith (P Int Rational)
sd' [IndexedPolynomialWith (P Int Rational)]
prs') [Int
1 ..] [P Int Rational]
qs

  -- Ignore log terms which end up being multiples of 0 = log 1.
  [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
-> Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return ([(P Int Rational, IndexedPolynomialWith (P Int Rational))]
 -> Maybe
      [(P Int Rational, IndexedPolynomialWith (P Int Rational))])
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
-> Maybe [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
forall a b. (a -> b) -> a -> b
$ ((P Int Rational, IndexedPolynomialWith (P Int Rational)) -> Bool)
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
-> [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
forall a. (a -> Bool) -> [a] -> [a]
filter (IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational) -> Bool
forall a. Eq a => a -> a -> Bool
(/=) IndexedPolynomialWith (P Int Rational)
1 (IndexedPolynomialWith (P Int Rational) -> Bool)
-> ((P Int Rational, IndexedPolynomialWith (P Int Rational))
    -> IndexedPolynomialWith (P Int Rational))
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (P Int Rational)
forall a b. (a, b) -> b
snd) [(P Int Rational, IndexedPolynomialWith (P Int Rational))]
terms
  where
    toTerm ::
      IndexedPolynomialWith IndexedPolynomial ->
      [IndexedPolynomialWith IndexedPolynomial] ->
      Int ->
      IndexedPolynomial ->
      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
    toTerm :: IndexedPolynomialWith (P Int Rational)
-> [IndexedPolynomialWith (P Int Rational)]
-> Int
-> P Int Rational
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
toTerm IndexedPolynomialWith (P Int Rational)
sd [IndexedPolynomialWith (P Int Rational)]
prs Int
i P Int Rational
q
      | P Int Rational -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree P Int Rational
q Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
1)
      | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== P Int Rational -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree P Int Rational
d = (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
sd)
      | (Just IndexedPolynomialWith (P Int Rational)
r) <- (IndexedPolynomialWith (P Int Rational) -> Bool)
-> [IndexedPolynomialWith (P Int Rational)]
-> Maybe (IndexedPolynomialWith (P Int Rational))
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
(==) Int
i (Int -> Bool)
-> (IndexedPolynomialWith (P Int Rational) -> Int)
-> IndexedPolynomialWith (P Int Rational)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith (P Int Rational) -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree) [IndexedPolynomialWith (P Int Rational)]
prs = P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
derive P Int Rational
q IndexedPolynomialWith (P Int Rational)
r
      | Bool
otherwise = (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
1)

    derive ::
      IndexedPolynomial ->
      IndexedPolynomialWith IndexedPolynomial ->
      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
    derive :: P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
derive P Int Rational
q IndexedPolynomialWith (P Int Rational)
s = (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
s')
      where
        as :: [P Int Rational]
as = P Int Rational -> [P Int Rational]
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 Int Rational -> [P Int Rational])
-> P Int Rational -> [P Int Rational]
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (P Int Rational) -> P Int Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient IndexedPolynomialWith (P Int Rational)
s
        s' :: IndexedPolynomialWith (P Int Rational)
s' = (IndexedPolynomialWith (P Int Rational)
 -> (Int, P Int Rational) -> IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (P Int Rational)
-> [(Int, P Int Rational)]
-> IndexedPolynomialWith (P Int Rational)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl IndexedPolynomialWith (P Int Rational)
-> (Int, P Int Rational) -> IndexedPolynomialWith (P Int Rational)
forall {p :: * -> * -> *} {e} {p :: * -> * -> *} {b}.
(Polynomial p e (P Int Rational), Polynomial p e (P Int Rational),
 Integral b, Num (p e (P Int Rational))) =>
p e (P Int Rational) -> (b, P Int Rational) -> p e (P Int Rational)
scalePoly IndexedPolynomialWith (P Int Rational)
s ([Int] -> [P Int Rational] -> [(Int, P Int Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int
1 ..] :: [Int]) [P Int Rational]
as)
          where
            scalePoly :: p e (P Int Rational) -> (b, P Int Rational) -> p e (P Int Rational)
scalePoly p e (P Int Rational)
x (b
j, P Int Rational
u) =
              Sum (p e (P Int Rational)) -> p e (P Int Rational)
forall a. Sum a -> a
getSum (Sum (p e (P Int Rational)) -> p e (P Int Rational))
-> Sum (p e (P Int Rational)) -> p e (P Int Rational)
forall a b. (a -> b) -> a -> b
$ (e -> P Int Rational -> Sum (p e (P Int Rational)))
-> p e (P Int Rational) -> Sum (p e (P Int Rational))
forall m.
Monoid m =>
(e -> P Int Rational -> m) -> p e (P Int Rational) -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (P Int Rational -> e -> P Int Rational -> Sum (p e (P Int Rational))
forall {p :: * -> * -> *} {e} {p :: * -> * -> *} {e} {c}.
(Polynomial p e (p e c), Polynomial p e c, Fractional c,
 Eq (p e c)) =>
p e c -> e -> p e c -> Sum (p e (p e c))
reduceTerm (P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (P Int Rational -> P Int Rational)
-> P Int Rational -> P Int Rational
forall a b. (a -> b) -> a -> b
$ P Int Rational -> P Int Rational -> P Int Rational
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 Int Rational
u P Int Rational
q P Int Rational -> b -> P Int Rational
forall a b. (Num a, Integral b) => a -> b -> a
^ b
j)) p e (P Int Rational)
x
            reduceTerm :: p e c -> e -> p e c -> Sum (p e (p e c))
reduceTerm p e c
v 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
$ p e c -> p e (p e c) -> p e (p e c)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (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)) =>
p e c -> p e c -> p e c
exactDivide p e c
c p e c
v) (p e (p e c) -> p e (p e c)) -> p e (p e c) -> p e (p e c)
forall a b. (a -> b) -> a -> b
$ e -> p e (p e c)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
e
            exactDivide :: p e c -> p e c -> p e c
exactDivide p e c
u p e c
v = p e c
r
              where
                (p e c
r, p e c
_) = p e c
u 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
v

-- | Given polynomials \(A\) and \(B\),
-- return a sum \(f\) of inverse tangents such that the following is true.
--
-- \[
-- \frac{df}{dx} = \frac{d}{dx} i \log \left( \frac{A + iB}{A - iB} \right)
-- \]
--
-- This allows integrals to be evaluated with only real-valued functions.
-- It also avoids the discontinuities in real-valued indefinite integrals which may result
-- when the integral uses logarithms with complex arguments.
--
-- For example,
--
-- >>> toHaskell $ simplify $ complexLogTermToAtan "x" (power 3 - 3 * power 1) (power 2 - 2)
-- "2 * atan x + 2 * atan ((x + (-3) * x ** 3 + x ** 5) / 2) + 2 * atan (x ** 3)"
--
-- so it is the case that
--
-- \[ \frac{d}{dx} \left( i \log \left( \frac{(x^3-3x) + i(x^2-2)}{(x^3-3x) - i(x^2-2)} \right) \right) =
-- \frac{d}{dx} \left( 2 \tan^{-1} \left(\frac{x^5-3x^3+x}{2}\right) + 2 \tan^{-1} \left(x^3\right) + 2 \tan^{-1} x \right) \]
complexLogTermToAtan ::
  -- | Symbol for the variable.
  Text ->
  -- | Polynomial \(A\).
  IndexedPolynomial ->
  -- | Polynomial \(B\).
  IndexedPolynomial ->
  -- | Sum \(f\) of inverse tangents.
  Expression
complexLogTermToAtan :: Text -> P Int Rational -> P Int Rational -> Expression
complexLogTermToAtan Text
v P Int Rational
a P Int Rational
b
  | P Int Rational
r P Int Rational -> P Int Rational -> Bool
forall a. Eq a => a -> a -> Bool
== P Int Rational
0 = Expression
2 Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression -> Expression
forall a. Floating a => a -> a
atan (Expression
a' Expression -> Expression -> Expression
forall a. Fractional a => a -> a -> a
/ Expression
b')
  | P Int Rational -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree P Int Rational
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< P Int Rational -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree P Int Rational
b = Text -> P Int Rational -> P Int Rational -> Expression
complexLogTermToAtan Text
v (-P Int Rational
b) P Int Rational
a
  | Bool
otherwise = Expression
2 Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression -> Expression
forall a. Floating a => a -> a
atan (Expression
s' Expression -> Expression -> Expression
forall a. Fractional a => a -> a -> a
/ Expression
g') Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
+ Text -> P Int Rational -> P Int Rational -> Expression
complexLogTermToAtan Text
v P Int Rational
d P Int Rational
c
  where
    (P Int Rational
_, P Int Rational
r) = P Int Rational
a P Int Rational
-> P Int Rational -> (P Int Rational, P Int Rational)
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 Int Rational
b
    (P Int Rational
d, P Int Rational
c, P Int Rational
g) = P Int Rational
-> P Int Rational
-> (P Int Rational, P Int Rational, P Int Rational)
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, p e c)
extendedEuclidean P Int Rational
b (-P Int Rational
a)
    a' :: Expression
a' = Text -> (Rational -> Expression) -> P Int Rational -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Rational -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient P Int Rational
a
    b' :: Expression
b' = Text -> (Rational -> Expression) -> P Int Rational -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Rational -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient P Int Rational
b
    g' :: Expression
g' = Text -> (Rational -> Expression) -> P Int Rational -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Rational -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient P Int Rational
g
    s' :: Expression
s' = Text -> (Rational -> Expression) -> P Int Rational -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Rational -> Expression
forall c. Real c => c -> Expression
toRationalCoefficient (P Int Rational -> Expression) -> P Int Rational -> Expression
forall a b. (a -> b) -> a -> b
$ P Int Rational
a P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
* P Int Rational
d P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
+ P Int Rational
b P Int Rational -> P Int Rational -> P Int Rational
forall a. Num a => a -> a -> a
* P Int Rational
c

-- | For the ingredients of a complex logarithm, return the ingredients of an equivalent real function in terms of an indefinite integral.
--
-- Specifically, for polynomials \(\left(R(t), S(t,x)\right)\) such that
--
-- \[
-- \frac{df}{dx} = \frac{d}{dx} \sum_{\alpha \in \{ t \mid R(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)
-- \]
--
-- then with return value \(\left( \left(P(u,v), Q(u,v)\right), \left(A(u,v,x), B(u,v,x)\right) \right)\),
-- and a return value \(g_{uv}\) from 'complexLogTermToAtan' for \(A(u,v)\) and \(B(u,v)\), the real function is
--
-- \[
-- \frac{df}{dx} = \frac{d}{dx} \left(
-- \sum_{(a,b) \in \{(u,v) \in (\mathbb{R}, \mathbb{R}) \mid P(u,v)=Q(u,v)=0, b > 0\}}
--   \left( a \log \left( A(a,b,x)^2 + B(a,b,x)^2 \right) + b g_{ab}(x) \right)
-- + \sum_{a \in \{t \in \mathbb{R} \mid R(t)=0 \}} \left( a \log (S(a,x)) \right)
-- \right)
-- \]
--
-- The return value are polynomials \(\left( (P,Q), (A,B) \right)\), where
--
-- * \(P\) is a \(u\)-polynomial, i.e., a polynomial with variable \(u\), with coefficients which are \(v\)-polynomials.
--
-- * \(Q\) is a \(u\)-polynomial, with coefficients which are \(v\)-polynomials.
--
-- * \(A\) is an \(x\)-polynomial, with coefficients which are \(u\)-polynomials, which in turn have coefficients with \(v\)-polynomials.
--
-- * \(B\) is an \(x\)-polynomial, with coefficients which are \(u\)-polynomials, which in turn have coefficients with \(v\)-polynomials.
--
-- For example,
--
-- >>> let r = 4 * power 2 + 1 :: IndexedPolynomial
-- >>> let s = power 3 + scale (2 * power 1) (power 2) - 3 * power 1 - scale (4 * power 1) 1 :: IndexedPolynomialWith IndexedPolynomial
-- >>> complexLogTermToRealTerm (r, s)
-- (([(0,(-4)x^2 + 1),(2,4)],[(1,8x)]),([(0,[(1,(-4))]),(1,[(0,(-3))]),(2,[(1,2)]),(3,[(0,1)])],[(0,[(0,(-4)x)]),(2,[(0,2x)])]))
--
-- While the return value may be hard to parse, this means:
--
-- \[
-- \begin{align*}
-- P & = 4u^2 - 4v^2 + 1 \\
-- Q & = 8uv \\
-- A & = x^3 + 2ux^2 - 3x - 4u \\
-- B & = 2vx^2 - 4v
-- \end{align*}
-- \]
complexLogTermToRealTerm ::
  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->
  ( (IndexedPolynomialWith IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial),
    (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial), IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
  )
complexLogTermToRealTerm :: (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> ((IndexedPolynomialWith (P Int Rational),
     IndexedPolynomialWith (P Int Rational)),
    (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)),
     IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
complexLogTermToRealTerm (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
s) = ((IndexedPolynomialWith (P Int Rational)
qp, IndexedPolynomialWith (P Int Rational)
qq), (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
sp, IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
sq))
  where
    -- For all of the following, i is the imaginary number.
    -- We use an i polynomial instead of Complex to represent complex numbers
    -- because the Complex a is not an instance of the Num class unless a is
    -- an instance of the RealFloat class.

    -- We use polynomial coefficients to introduce a separate variable.
    -- An alternative would have been to use Expression coefficients,
    -- but this would require a guarantee that we can rewrite an Expression
    -- down to the degree where we can tease apart the real and imaginary parts
    -- in a complex number.

    -- Compute q(u+iv) as an i polynomial with coefficients
    -- of u polynomials with coefficients
    -- of v polynomials with rational coefficients.
    q' :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
q' = Sum
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
 -> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomialWith (P Int Rational)
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall m.
Monoid m =>
(Int -> IndexedPolynomialWith (P Int Rational) -> m)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomialWith (P Int Rational)
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (a -> b) -> a -> b
$ Sum
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
 -> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> P Int Rational
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall m. Monoid m => (Int -> Rational -> m) -> P Int Rational -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> Rational
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
fromTerm P Int Rational
q
      where
        fromTerm :: Int -> Rational -> Sum (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
        fromTerm :: Int
-> Rational
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
fromTerm Int
e Rational
c = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. a -> Sum a
Sum (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
c' IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
u IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
i IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
v) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
c' = IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational -> P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c P Int Rational
1) IndexedPolynomialWith (P Int Rational)
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1
        i :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
i = Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        u :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
u = IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1
        v :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
v = IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> P Int Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (P Int Rational)
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1
    -- q' == qp + i * qq
    (IndexedPolynomialWith (P Int Rational)
qp, IndexedPolynomialWith (P Int Rational)
qq) = (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Int -> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
q' Int
0, IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Int -> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
q' Int
1)

    -- Compute s(u+iv,x) as an i polynomial with coefficients
    -- of x polynomials with coefficients
    -- of u polynomials with coefficients
    -- of v polynomials with rational coefficients.
    s' :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
s' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall m.
Monoid m =>
(Int
 -> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
 -> m)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a b. (a -> b) -> a -> b
$ Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (a -> b) -> a -> b
$ (Int
 -> P Int Rational
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> IndexedPolynomialWith (P Int Rational)
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall m.
Monoid m =>
(Int -> P Int Rational -> m)
-> IndexedPolynomialWith (P Int Rational) -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> P Int Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
fromTerm IndexedPolynomialWith (P Int Rational)
s
      where
        fromTerm :: Int -> IndexedPolynomial -> Sum (IndexedPolynomialWith (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
        fromTerm :: Int
-> P Int Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
fromTerm Int
e P Int Rational
c = IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
c' IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
x IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
c' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> P Int Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall m. Monoid m => (Int -> Rational -> m) -> P Int Rational -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall {b}.
Integral b =>
b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
fromCoefficient P Int Rational
c
            fromCoefficient :: b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
fromCoefficient b
e' Rational
c'' = IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
c''' IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
u IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
i IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
v) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> b
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall a b. (Num a, Integral b) => a -> b -> a
^ b
e'
              where
                c''' :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
c''' = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational -> P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c'' P Int Rational
1) IndexedPolynomialWith (P Int Rational)
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
1
        i :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
i = Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        x :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
x = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
1
        u :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
u = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
1
        v :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
v = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int Rational
-> IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> P Int Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (P Int Rational)
1) IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
1
    -- s' = sp + i * sq
    (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
sp, IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
sq) = (IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
s' Int
0, IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)))
s' Int
1)

    -- For terms in polynomials of i, reduce them to the form x or i*x.
    reduceImaginary :: (Eq a, Num a) => Int -> a -> Sum (IndexedPolynomialWith a)
    reduceImaginary :: forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary Int
e a
c = IndexedPolynomialWith a -> Sum (IndexedPolynomialWith a)
forall a. a -> Sum a
Sum (IndexedPolynomialWith a -> Sum (IndexedPolynomialWith a))
-> IndexedPolynomialWith a -> Sum (IndexedPolynomialWith a)
forall a b. (a -> b) -> a -> b
$ case Int
e Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
4 of
      Int
0 -> IndexedPolynomialWith a
c'
      Int
1 -> IndexedPolynomialWith a
c' IndexedPolynomialWith a
-> IndexedPolynomialWith a -> IndexedPolynomialWith a
forall a. Num a => a -> a -> a
* IndexedPolynomialWith a
i
      Int
2 -> IndexedPolynomialWith a
c' IndexedPolynomialWith a
-> IndexedPolynomialWith a -> IndexedPolynomialWith a
forall a. Num a => a -> a -> a
* (-IndexedPolynomialWith a
1)
      Int
3 -> IndexedPolynomialWith a
c' IndexedPolynomialWith a
-> IndexedPolynomialWith a -> IndexedPolynomialWith a
forall a. Num a => a -> a -> a
* (-IndexedPolynomialWith a
i)
      Int
_ -> IndexedPolynomialWith a
0 -- Not possible.
      where
        i :: IndexedPolynomialWith a
i = Int -> IndexedPolynomialWith a
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        c' :: IndexedPolynomialWith a
c' = a -> IndexedPolynomialWith a -> IndexedPolynomialWith a
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale a
c IndexedPolynomialWith a
1

-- | For the ingredients of a complex logarithm, return an equivalent real function in terms of an indefinite integral.
--
-- Specifically, for polynomials \(\left(R(t), S(t,x)\right)\) such that
--
-- \[
-- \frac{df}{dx} = \frac{d}{dx} \sum_{\alpha \in \{ t \mid R(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)
-- \]
--
-- a symbolic representation for \(f\) will be returned.  See 'complexLogTermToRealTerm' for specifics as to how \(f\) is derived.
complexLogTermToRealExpression ::
  -- | Symbol for the variable.
  Text ->
  -- | Polynomials \(R(t)\) and \(S(t,x)\).
  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->
  -- | Expression for the real function \(f\).
  Maybe Expression
complexLogTermToRealExpression :: Text
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> Maybe Expression
complexLogTermToRealExpression Text
v (P Int Rational
r, IndexedPolynomialWith (P Int Rational)
s)
  | (Just [(Rational, Rational)]
xys) <- IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
-> Maybe [(Rational, Rational)]
solveBivariatePolynomials IndexedPolynomialWith (P Int Rational)
p IndexedPolynomialWith (P Int Rational)
q,
    (Just [Expression]
h) <- [(Rational, Rational)] -> Maybe [Expression]
f [(Rational, Rational)]
xys,
    (Just [Rational]
zs) <- Maybe [Expression] -> Maybe [Rational]
toRationalList (P Int Rational -> Maybe [Expression]
solve P Int Rational
r) =
      Expression -> Maybe Expression
forall a. a -> Maybe a
Just (Expression -> Maybe Expression) -> Expression -> Maybe Expression
forall a b. (a -> b) -> a -> b
$ [Expression] -> Expression
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Expression]
h Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
+ [Rational] -> Expression
forall {t :: * -> *}.
(Foldable t, Monad t) =>
t Rational -> Expression
g [Rational]
zs
  | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing
  where
    ((IndexedPolynomialWith (P Int Rational)
p, IndexedPolynomialWith (P Int Rational)
q), (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
a, IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
b)) = (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> ((IndexedPolynomialWith (P Int Rational),
     IndexedPolynomialWith (P Int Rational)),
    (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational)),
     IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))))
complexLogTermToRealTerm (P Int Rational
r, IndexedPolynomialWith (P Int Rational)
s)

    f :: [(Rational, Rational)] -> Maybe [Expression]
    f :: [(Rational, Rational)] -> Maybe [Expression]
f [(Rational, Rational)]
xys = [Maybe Expression] -> Maybe [Expression]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence ([Maybe Expression] -> Maybe [Expression])
-> [Maybe Expression] -> Maybe [Expression]
forall a b. (a -> b) -> a -> b
$ do
      (Rational
x, Rational
y) <- ((Rational, Rational) -> Bool)
-> [(Rational, Rational)] -> [(Rational, Rational)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
0) (Rational -> Bool)
-> ((Rational, Rational) -> Rational)
-> (Rational, Rational)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
forall a b. (a, b) -> b
snd) [(Rational, Rational)]
xys
      let flatten'' :: IndexedPolynomialWith (P Int Rational) -> P Int Expression
flatten'' = (P Int Rational -> Expression)
-> IndexedPolynomialWith (P Int Rational) -> P Int Expression
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 (Expression
-> (Rational -> Expression) -> P Int Rational -> Expression
forall {p :: * -> * -> *} {a} {t}.
Polynomial p a t =>
Expression -> (t -> Expression) -> p a t -> Expression
toExpr (Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
y) Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational) -- v-polynomials into Expressions.
      let flatten' :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> P Int Expression
flatten' = (IndexedPolynomialWith (P Int Rational) -> Expression)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> P Int Expression
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 (Expression
-> (Expression -> Expression) -> P Int Expression -> Expression
forall {p :: * -> * -> *} {a} {t}.
Polynomial p a t =>
Expression -> (t -> Expression) -> p a t -> Expression
toExpr (Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
x) Expression -> Expression
forall a. a -> a
id (P Int Expression -> Expression)
-> (IndexedPolynomialWith (P Int Rational) -> P Int Expression)
-> IndexedPolynomialWith (P Int Rational)
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith (P Int Rational) -> P Int Expression
flatten'') -- u-polynomials into Expressions.
      let flatten :: IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Expression
flatten = Expression
-> (Expression -> Expression) -> P Int Expression -> Expression
forall {p :: * -> * -> *} {a} {t}.
Polynomial p a t =>
Expression -> (t -> Expression) -> p a t -> Expression
toExpr (Text -> Expression
Symbol Text
v) Expression -> Expression
forall a. a -> a
id (P Int Expression -> Expression)
-> (IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
    -> P Int Expression)
-> IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> P Int Expression
flatten' -- x-polynomials into Expressions.
      -- a and b flattened into Expressions.
      let a' :: Expression
a' = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Expression
flatten IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
a
      let b' :: Expression
b' = IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> Expression
flatten IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
b
      -- a and b flattened into x-polynomials with rational number coefficients.
      Maybe Expression -> [Maybe Expression]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe Expression -> [Maybe Expression])
-> Maybe Expression -> [Maybe Expression]
forall a b. (a -> b) -> a -> b
$ do
        P Int Rational
a'' <- P Int Expression -> Maybe (P Int Rational)
convertCoefficients (P Int Expression -> Maybe (P Int Rational))
-> P Int Expression -> Maybe (P Int Rational)
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> P Int Expression
flatten' IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
a
        P Int Rational
b'' <- P Int Expression -> Maybe (P Int Rational)
convertCoefficients (P Int Expression -> Maybe (P Int Rational))
-> P Int Expression -> Maybe (P Int Rational)
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
-> P Int Expression
flatten' IndexedPolynomialWith (IndexedPolynomialWith (P Int Rational))
b
        Expression -> Maybe Expression
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (Expression -> Maybe Expression) -> Expression -> Maybe Expression
forall a b. (a -> b) -> a -> b
$ Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
x Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression -> Expression
forall a. Floating a => a -> a
log (Expression
a' Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression
a' Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
+ Expression
b' Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression
b') Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
+ Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
y Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Text -> P Int Rational -> P Int Rational -> Expression
complexLogTermToAtan Text
v P Int Rational
a'' P Int Rational
b''

    g :: t Rational -> Expression
g t Rational
zs = t Expression -> Expression
forall a. Num a => t a -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (t Expression -> Expression) -> t Expression -> Expression
forall a b. (a -> b) -> a -> b
$ do
      Rational
z <- t Rational
zs
      let s' :: P Int Expression
s' = (P Int Rational -> Expression)
-> IndexedPolynomialWith (P Int Rational) -> P Int Expression
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 (Expression
-> (Rational -> Expression) -> P Int Rational -> Expression
forall {p :: * -> * -> *} {a} {t}.
Polynomial p a t =>
Expression -> (t -> Expression) -> p a t -> Expression
toExpr (Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
z) Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational) IndexedPolynomialWith (P Int Rational)
s
      Expression -> t Expression
forall a. a -> t a
forall (m :: * -> *) a. Monad m => a -> m a
return (Expression -> t Expression) -> Expression -> t Expression
forall a b. (a -> b) -> a -> b
$ Rational -> Expression
forall a. Fractional a => Rational -> a
fromRational Rational
z Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression -> Expression
Log' (Text
-> (Expression -> Expression) -> P Int Expression -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Expression -> Expression
toSymbolicCoefficient P Int Expression
s')

    toRationalList :: Maybe [Expression] -> Maybe [Rational]
    toRationalList :: Maybe [Expression] -> Maybe [Rational]
toRationalList Maybe [Expression]
Nothing = Maybe [Rational]
forall a. Maybe a
Nothing
    toRationalList (Just []) = [Rational] -> Maybe [Rational]
forall a. a -> Maybe a
Just []
    toRationalList (Just (Expression
x : [Expression]
xs))
      | (Just Rational
x'') <- Expression -> Maybe Rational
forall {a}. Fractional a => Expression -> Maybe a
convert (Expression -> Expression
simplify Expression
x'), (Just [Rational]
xs'') <- Maybe [Rational]
xs' = [Rational] -> Maybe [Rational]
forall a. a -> Maybe a
Just ([Rational] -> Maybe [Rational]) -> [Rational] -> Maybe [Rational]
forall a b. (a -> b) -> a -> b
$ Rational
x'' Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [Rational]
xs''
      | Bool
otherwise = Maybe [Rational]
forall a. Maybe a
Nothing
      where
        x' :: Expression
x' = Expression -> Expression
simplify Expression
x
        xs' :: Maybe [Rational]
xs' = Maybe [Expression] -> Maybe [Rational]
toRationalList (Maybe [Expression] -> Maybe [Rational])
-> Maybe [Expression] -> Maybe [Rational]
forall a b. (a -> b) -> a -> b
$ [Expression] -> Maybe [Expression]
forall a. a -> Maybe a
Just [Expression]
xs

    -- Convert a simplified Expression into a rational number.
    convert :: Expression -> Maybe a
convert (Number Integer
n) = a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n
    convert (Number Integer
n :/: Number Integer
m) = a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n a -> a -> a
forall a. Fractional a => a -> a -> a
/ Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
m
    convert Expression
_ = Maybe a
forall a. Maybe a
Nothing

    -- Convert polynomial with Expression coefficients into a polynomial with rational number coefficients.
    convertCoefficients :: IndexedPolynomialWith Expression -> Maybe IndexedPolynomial
    convertCoefficients :: P Int Expression -> Maybe (P Int Rational)
convertCoefficients P Int Expression
x = [P Int Rational] -> P Int Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([P Int Rational] -> P Int Rational)
-> ([(Int, Rational)] -> [P Int Rational])
-> [(Int, Rational)]
-> P Int Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Rational) -> P Int Rational)
-> [(Int, Rational)] -> [P Int Rational]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
e, Rational
c) -> Rational -> P Int Rational -> P Int Rational
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c (Int -> P Int Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
e)) ([(Int, Rational)] -> P Int Rational)
-> Maybe [(Int, Rational)] -> Maybe (P Int Rational)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Maybe (Int, Rational)] -> Maybe [(Int, Rational)]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence ((Int -> Expression -> [Maybe (Int, Rational)])
-> P Int Expression -> [Maybe (Int, Rational)]
forall m.
Monoid m =>
(Int -> Expression -> m) -> P Int Expression -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\Int
e Expression
c -> [(Int
e,) (Rational -> (Int, Rational))
-> Maybe Rational -> Maybe (Int, Rational)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expression -> Maybe Rational
forall {a}. Fractional a => Expression -> Maybe a
convert (Expression -> Expression
simplify Expression
c)]) P Int Expression
x)

    -- Turns a polynomial into an Expression.
    -- Function h is used to turn the coefficient into an Expression.
    toExpr :: Expression -> (t -> Expression) -> p a t -> Expression
toExpr Expression
x t -> Expression
h p a t
u = Sum Expression -> Expression
forall a. Sum a -> a
getSum (Sum Expression -> Expression) -> Sum Expression -> Expression
forall a b. (a -> b) -> a -> b
$ (a -> t -> Sum Expression) -> p a t -> Sum Expression
forall m. Monoid m => (a -> t -> m) -> p a t -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\a
e'' t
c -> Expression -> Sum Expression
forall a. a -> Sum a
Sum (Expression -> Sum Expression) -> Expression -> Sum Expression
forall a b. (a -> b) -> a -> b
$ t -> Expression
h t
c Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* (Expression
x Expression -> Expression -> Expression
forall a. Floating a => a -> a -> a
** Integer -> Expression
Number (a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
e''))) p a t
u

-- | From the ingredients of a complex logarithm, return the expression for the complex algorithm.
-- Specifically, for polynomials \(\left(Q(t), S(t,x)\right)\),
-- a symbolic representation for the following will be returned.
--
-- \[
-- \sum_{\alpha \in \{ t \mid Q(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)
-- \]
complexLogTermToComplexExpression ::
  -- | Symbol for the variable.
  Text ->
  -- | Polynomials \(Q(t)\) and \(S(t,x)\).
  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->
  -- | Expression for the logarithm.
  Maybe Expression
complexLogTermToComplexExpression :: Text
-> (P Int Rational, IndexedPolynomialWith (P Int Rational))
-> Maybe Expression
complexLogTermToComplexExpression Text
v (P Int Rational
q, IndexedPolynomialWith (P Int Rational)
s) = do
  [Expression]
as <- P Int Rational -> Maybe [Expression]
complexSolve P Int Rational
q
  let terms :: [Expression]
terms = do
        Expression
a <- [Expression]
as
        let s' :: P Int Expression
s' = (P Int Rational -> Expression)
-> IndexedPolynomialWith (P Int Rational) -> P Int Expression
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 (Expression -> P Int Rational -> Expression
forall {p :: * -> * -> *} {a} {a}.
(Polynomial p a Rational, Floating a) =>
a -> p a Rational -> a
collapse Expression
a) IndexedPolynomialWith (P Int Rational)
s
        let s'' :: Expression
s'' = Text
-> (Expression -> Expression) -> P Int Expression -> Expression
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
Text -> (c -> Expression) -> p e c -> Expression
toExpression Text
v Expression -> Expression
toSymbolicCoefficient P Int Expression
s'
        Expression -> [Expression]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return (Expression -> [Expression]) -> Expression -> [Expression]
forall a b. (a -> b) -> a -> b
$ Expression
a Expression -> Expression -> Expression
forall a. Num a => a -> a -> a
* Expression -> Expression
forall a. Floating a => a -> a
log Expression
s''
  Expression -> Maybe Expression
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (Expression -> Maybe Expression) -> Expression -> Maybe Expression
forall a b. (a -> b) -> a -> b
$ [Expression] -> Expression
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Expression]
terms
  where
    -- Collapse a polynomial coefficient of a polynomial into an expression with the variable substituted.
    -- E.g., turn (t+2)x+1 into (3+2)x+1 for t=3.
    collapse :: a -> p a Rational -> a
collapse a
a p a Rational
c' = Sum a -> a
forall a. Sum a -> a
getSum (Sum a -> a) -> Sum a -> a
forall a b. (a -> b) -> a -> b
$ (a -> Rational -> Sum a) -> p a Rational -> Sum a
forall m. Monoid m => (a -> Rational -> m) -> p a Rational -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\a
e Rational
c -> a -> Sum a
forall a. a -> Sum a
Sum (a -> Sum a) -> a -> Sum a
forall a b. (a -> b) -> a -> b
$ Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
c a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Floating a => a -> a -> a
** a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
e) p a Rational
c'

-- | Returns the roots for two variables in two polynomials.
--
-- Only supports rational roots.  If not all real roots are rational, then it will return 'Nothing'.
-- Returning all real roots would be preferable, but this is not supported at this time.
--
-- If the function cannot derive the roots otherwise, either, 'Nothing' will be returned as well.
solveBivariatePolynomials ::
  IndexedPolynomialWith IndexedPolynomial ->
  IndexedPolynomialWith IndexedPolynomial ->
  Maybe [(Rational, Rational)]
solveBivariatePolynomials :: IndexedPolynomialWith (P Int Rational)
-> IndexedPolynomialWith (P Int Rational)
-> Maybe [(Rational, Rational)]
solveBivariatePolynomials IndexedPolynomialWith (P Int Rational)
p IndexedPolynomialWith (P Int Rational)
q = do
  let p' :: P Int (Function (P Int Rational))
p' = IndexedPolynomialWith (P Int Rational)
-> P Int (Function (P Int Rational))
toRationalFunctionCoefficients IndexedPolynomialWith (P Int Rational)
p
  let q' :: P Int (Function (P Int Rational))
q' = IndexedPolynomialWith (P Int Rational)
-> P Int (Function (P Int Rational))
toRationalFunctionCoefficients IndexedPolynomialWith (P Int Rational)
q
  P Int Rational
resultant <- Function (P Int Rational) -> Maybe (P Int Rational)
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 Int Rational) -> Maybe (P Int Rational))
-> Function (P Int Rational) -> Maybe (P Int Rational)
forall a b. (a -> b) -> a -> b
$ (Function (P Int Rational), [P Int (Function (P Int Rational))])
-> Function (P Int Rational)
forall a b. (a, b) -> a
fst ((Function (P Int Rational), [P Int (Function (P Int Rational))])
 -> Function (P Int Rational))
-> (Function (P Int Rational), [P Int (Function (P Int Rational))])
-> Function (P Int Rational)
forall a b. (a -> b) -> a -> b
$ P Int (Function (P Int Rational))
-> P Int (Function (P Int Rational))
-> (Function (P Int Rational), [P Int (Function (P Int Rational))])
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 Int (Function (P Int Rational))
p' P Int (Function (P Int Rational))
q'
  [Expression]
vs' <- P Int Rational -> Maybe [Expression]
solve P Int Rational
resultant
  [Rational]
vs <- (Expression -> Maybe Rational) -> [Expression] -> Maybe [Rational]
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 (Expression -> Maybe Rational
convert (Expression -> Maybe Rational)
-> (Expression -> Expression) -> Expression -> Maybe Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression
simplify) [Expression]
vs'
  [[(Rational, Rational)]] -> [(Rational, Rational)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[(Rational, Rational)]] -> [(Rational, Rational)])
-> Maybe [[(Rational, Rational)]] -> Maybe [(Rational, Rational)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Rational -> Maybe [(Rational, Rational)])
-> [Rational] -> Maybe [[(Rational, Rational)]]
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 Rational -> Maybe [(Rational, Rational)]
solveForU [Rational]
vs
  where
    toRationalFunctionCoefficients :: IndexedPolynomialWith (P Int Rational)
-> P Int (Function (P Int Rational))
toRationalFunctionCoefficients = (P Int Rational -> Function (P Int Rational))
-> IndexedPolynomialWith (P Int Rational)
-> P Int (Function (P Int Rational))
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 Int Rational -> Function (P Int Rational)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial

    -- For each v, returns list of (u,v) such that P(u,v)=Q(u,v)=0.
    solveForU :: Rational -> Maybe [(Rational, Rational)]
    solveForU :: Rational -> Maybe [(Rational, Rational)]
solveForU Rational
v
      | P Int Rational
0 <- P Int Rational
p' = do
          -- Any u will make p'=0 true, so we only need to solve p'.
          [Maybe Rational]
u <- (Expression -> Maybe Rational) -> [Expression] -> [Maybe Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Expression -> Maybe Rational
convert (Expression -> Maybe Rational)
-> (Expression -> Expression) -> Expression -> Maybe Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression
simplify) ([Expression] -> [Maybe Rational])
-> Maybe [Expression] -> Maybe [Maybe Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P Int Rational -> Maybe [Expression]
solve P Int Rational
q'
          (Rational -> (Rational, Rational))
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map (,Rational
v) ([Rational] -> [(Rational, Rational)])
-> Maybe [Rational] -> Maybe [(Rational, Rational)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Maybe Rational] -> Maybe [Rational]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence [Maybe Rational]
u
      | P Int Rational
0 <- P Int Rational
q' = do
          -- Any u will make q'=0 true, so we only need to solve p'.
          [Maybe Rational]
u <- (Expression -> Maybe Rational) -> [Expression] -> [Maybe Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Expression -> Maybe Rational
convert (Expression -> Maybe Rational)
-> (Expression -> Expression) -> Expression -> Maybe Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression
simplify) ([Expression] -> [Maybe Rational])
-> Maybe [Expression] -> Maybe [Maybe Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P Int Rational -> Maybe [Expression]
solve P Int Rational
p'
          (Rational -> (Rational, Rational))
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map (,Rational
v) ([Rational] -> [(Rational, Rational)])
-> Maybe [Rational] -> Maybe [(Rational, Rational)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Maybe Rational] -> Maybe [Rational]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence [Maybe Rational]
u
      | Bool
otherwise = do
          [Maybe Rational]
up <- (Expression -> Maybe Rational) -> [Expression] -> [Maybe Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Expression -> Maybe Rational
convert (Expression -> Maybe Rational)
-> (Expression -> Expression) -> Expression -> Maybe Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression
simplify) ([Expression] -> [Maybe Rational])
-> Maybe [Expression] -> Maybe [Maybe Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P Int Rational -> Maybe [Expression]
solve P Int Rational
p'
          [Maybe Rational]
uq <- (Expression -> Maybe Rational) -> [Expression] -> [Maybe Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Expression -> Maybe Rational
convert (Expression -> Maybe Rational)
-> (Expression -> Expression) -> Expression -> Maybe Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expression -> Expression
simplify) ([Expression] -> [Maybe Rational])
-> Maybe [Expression] -> Maybe [Maybe Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P Int Rational -> Maybe [Expression]
solve P Int Rational
q'
          [Rational]
up' <- [Maybe Rational] -> Maybe [Rational]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence [Maybe Rational]
up
          [Rational]
uq' <- [Maybe Rational] -> Maybe [Rational]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence [Maybe Rational]
uq
          [(Rational, Rational)] -> Maybe [(Rational, Rational)]
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return ([(Rational, Rational)] -> Maybe [(Rational, Rational)])
-> [(Rational, Rational)] -> Maybe [(Rational, Rational)]
forall a b. (a -> b) -> a -> b
$ (Rational -> (Rational, Rational))
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map (,Rational
v) ([Rational] -> [(Rational, Rational)])
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> a -> b
$ [Rational]
up' [Rational] -> [Rational] -> [Rational]
forall a. Eq a => [a] -> [a] -> [a]
`intersect` [Rational]
uq'
      where
        p' :: P Int Rational
p' = (P Int Rational -> Rational)
-> IndexedPolynomialWith (P Int Rational) -> P Int Rational
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 (Sum Rational -> Rational
forall a. Sum a -> a
getSum (Sum Rational -> Rational)
-> (P Int Rational -> Sum Rational) -> P Int Rational -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Rational -> Sum Rational) -> P Int Rational -> Sum Rational
forall m. Monoid m => (Int -> Rational -> m) -> P Int Rational -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\Int
e Rational
c -> Rational -> Sum Rational
forall a. a -> Sum a
Sum (Rational -> Sum Rational) -> Rational -> Sum Rational
forall a b. (a -> b) -> a -> b
$ Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
v Rational -> Int -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e)) IndexedPolynomialWith (P Int Rational)
p
        q' :: P Int Rational
q' = (P Int Rational -> Rational)
-> IndexedPolynomialWith (P Int Rational) -> P Int Rational
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 (Sum Rational -> Rational
forall a. Sum a -> a
getSum (Sum Rational -> Rational)
-> (P Int Rational -> Sum Rational) -> P Int Rational -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Rational -> Sum Rational) -> P Int Rational -> Sum Rational
forall m. Monoid m => (Int -> Rational -> m) -> P Int Rational -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\Int
e Rational
c -> Rational -> Sum Rational
forall a. a -> Sum a
Sum (Rational -> Sum Rational) -> Rational -> Sum Rational
forall a b. (a -> b) -> a -> b
$ Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
v Rational -> Int -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e)) IndexedPolynomialWith (P Int Rational)
q

    -- Turn a simplified Expression into a rational number if possible.
    convert :: Expression -> Maybe Rational
    convert :: Expression -> Maybe Rational
convert (Number Integer
n) = Rational -> Maybe Rational
forall a. a -> Maybe a
Just (Rational -> Maybe Rational) -> Rational -> Maybe Rational
forall a b. (a -> b) -> a -> b
$ Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n
    convert (Number Integer
n :/: Number Integer
m) = Rational -> Maybe Rational
forall a. a -> Maybe a
Just (Rational -> Maybe Rational) -> Rational -> Maybe Rational
forall a b. (a -> b) -> a -> b
$ Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
m
    convert Expression
_ = Maybe Rational
forall a. Maybe a
Nothing