-- |
-- 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,

    -- * Support

    -- | Functions and types useful when integrating rational functions.
    toRationalFunction,
    RationalFunction (..),
  )
where

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.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.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)\]
integrate :: Text -> Expression -> Maybe Expression
integrate :: Text -> Expression -> Maybe Expression
integrate Text
v Expression
e
  | (Expression
x :/: Expression
y) <- Expression
e',
    (Just IndexedPolynomial
n) <- (Text -> Maybe IndexedPolynomial, Expression -> Maybe Rational)
-> Expression -> Maybe IndexedPolynomial
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 IndexedPolynomial, 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 IndexedPolynomial
d) <- (Text -> Maybe IndexedPolynomial, Expression -> Maybe Rational)
-> Expression -> Maybe IndexedPolynomial
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 IndexedPolynomial, 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,
    IndexedPolynomial
d IndexedPolynomial -> IndexedPolynomial -> Bool
forall a. Eq a => a -> a -> Bool
/= IndexedPolynomial
0 =
      IndexedPolynomial -> IndexedPolynomial -> Maybe Expression
integrate' IndexedPolynomial
n IndexedPolynomial
d
  | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing
  where
    e' :: Expression
e' = Text -> Expression -> Expression
simplifyForVariable Text
v Expression
e
    integrate' :: IndexedPolynomial -> IndexedPolynomial -> Maybe Expression
integrate' IndexedPolynomial
n IndexedPolynomial
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.
        ([RationalFunction]
g, RationalFunction
h) = RationalFunction -> ([RationalFunction], RationalFunction)
hermiteReduce (RationalFunction -> ([RationalFunction], RationalFunction))
-> RationalFunction -> ([RationalFunction], RationalFunction)
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
n IndexedPolynomial
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
$ (RationalFunction -> Expression)
-> [RationalFunction] -> [Expression]
forall a b. (a -> b) -> [a] -> [b]
map RationalFunction -> Expression
fromRationalFunction [RationalFunction]
g

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

        -- Derive the log terms in the integral.
        h' :: RationalFunction
h' = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
r IndexedPolynomial
denom
        logTerms :: Maybe [(IndexedPolynomial, P Int IndexedPolynomial)]
logTerms = RationalFunction
-> Maybe [(IndexedPolynomial, P Int IndexedPolynomial)]
rationalIntegralLogTerms RationalFunction
h'
        logs :: Maybe Expression
        logs :: Maybe Expression
logs
          | (Just [(IndexedPolynomial, P Int IndexedPolynomial)]
terms) <- Maybe [(IndexedPolynomial, P Int IndexedPolynomial)]
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
<$> [Maybe Expression] -> Maybe [Expression]
forall a. [Maybe a] -> Maybe [a]
toMaybeList (((IndexedPolynomial, P Int IndexedPolynomial) -> Maybe Expression)
-> [(IndexedPolynomial, P Int IndexedPolynomial)]
-> [Maybe Expression]
forall a b. (a -> b) -> [a] -> [b]
map (Text
-> (IndexedPolynomial, P Int IndexedPolynomial) -> Maybe Expression
complexLogTermToRealExpression Text
v) [(IndexedPolynomial, P Int IndexedPolynomial)]
terms)
          | Bool
otherwise = Maybe Expression
forall a. Maybe a
Nothing

        fromRationalFunction :: RationalFunction -> Expression
fromRationalFunction (RationalFunction IndexedPolynomial
u IndexedPolynomial
w) = Expression
u' Expression -> Expression -> Expression
forall a. Fractional a => a -> a -> a
/ Expression
w'
          where
            u' :: Expression
u' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 IndexedPolynomial
u
            w' :: Expression
w' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 IndexedPolynomial
w

-- | Represents the ratio of two polynomials with rational number coefficients.
data RationalFunction = RationalFunction IndexedPolynomial IndexedPolynomial
  deriving (RationalFunction -> RationalFunction -> Bool
(RationalFunction -> RationalFunction -> Bool)
-> (RationalFunction -> RationalFunction -> Bool)
-> Eq RationalFunction
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RationalFunction -> RationalFunction -> Bool
== :: RationalFunction -> RationalFunction -> Bool
$c/= :: RationalFunction -> RationalFunction -> Bool
/= :: RationalFunction -> RationalFunction -> Bool
Eq)

instance Show RationalFunction where
  show :: RationalFunction -> String
show (RationalFunction IndexedPolynomial
n IndexedPolynomial
d) = String
"(" String -> ShowS
forall a. Semigroup a => a -> a -> a
<> IndexedPolynomial -> String
forall a. Show a => a -> String
show IndexedPolynomial
n String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
") / (" String -> ShowS
forall a. Semigroup a => a -> a -> a
<> IndexedPolynomial -> String
forall a. Show a => a -> String
show IndexedPolynomial
d String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
")"

-- | The numerator and denominator in the results
-- for '(+)', '(-)', '(*)', and 'negate' will be coprime.
instance Num RationalFunction where
  (RationalFunction IndexedPolynomial
x IndexedPolynomial
y) + :: RationalFunction -> RationalFunction -> RationalFunction
+ (RationalFunction IndexedPolynomial
u IndexedPolynomial
v) =
    IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction (IndexedPolynomial
x IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
v IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
+ IndexedPolynomial
u IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
y) (IndexedPolynomial
y IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
v)

  (RationalFunction IndexedPolynomial
x IndexedPolynomial
y) - :: RationalFunction -> RationalFunction -> RationalFunction
- (RationalFunction IndexedPolynomial
u IndexedPolynomial
v) =
    IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction (IndexedPolynomial
x IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
v IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
- IndexedPolynomial
u IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
y) (IndexedPolynomial
y IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
v)

  (RationalFunction IndexedPolynomial
x IndexedPolynomial
y) * :: RationalFunction -> RationalFunction -> RationalFunction
* (RationalFunction IndexedPolynomial
u IndexedPolynomial
v) =
    IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction (IndexedPolynomial
x IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
u) (IndexedPolynomial
y IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
v)

  abs :: RationalFunction -> RationalFunction
abs = RationalFunction -> RationalFunction
forall a. a -> a
id

  signum :: RationalFunction -> RationalFunction
signum RationalFunction
0 = RationalFunction
0
  signum RationalFunction
_ = RationalFunction
1

  fromInteger :: Integer -> RationalFunction
fromInteger Integer
n = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction (Integer -> IndexedPolynomial
forall a. Num a => Integer -> a
fromInteger Integer
n) IndexedPolynomial
1

instance Fractional RationalFunction where
  fromRational :: Rational -> RationalFunction
fromRational Rational
q = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction (Rational -> IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
q IndexedPolynomial
1) IndexedPolynomial
1
  recip :: RationalFunction -> RationalFunction
recip (RationalFunction IndexedPolynomial
p IndexedPolynomial
q) = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction IndexedPolynomial
q IndexedPolynomial
p

-- | Form a rational function from two polynomials.
-- The polynomials will be reduced so that the numerator and denominator are coprime.
toRationalFunction ::
  -- | Numerator.
  IndexedPolynomial ->
  -- | Denominator.
  IndexedPolynomial ->
  RationalFunction
toRationalFunction :: IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
x IndexedPolynomial
0 = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction IndexedPolynomial
x IndexedPolynomial
0
toRationalFunction IndexedPolynomial
x IndexedPolynomial
y = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction IndexedPolynomial
x' IndexedPolynomial
y'
  where
    g :: IndexedPolynomial
g = IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
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 IndexedPolynomial
x IndexedPolynomial
y
    (IndexedPolynomial
x', IndexedPolynomial
_) = IndexedPolynomial
x IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
g
    (IndexedPolynomial
y', IndexedPolynomial
_) = IndexedPolynomial
y IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
g

-- | 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 $ toRationalFunction p q
-- ([(3) / (x^2 + 2),(8x^2 + 4) / (x^5 + 4x^3 + 4x)],(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 :: RationalFunction -> ([RationalFunction], RationalFunction)
hermiteReduce h :: RationalFunction
h@(RationalFunction IndexedPolynomial
_ IndexedPolynomial
0) = ([], RationalFunction
h)
hermiteReduce h :: RationalFunction
h@(RationalFunction IndexedPolynomial
x IndexedPolynomial
y)
  | (Just ([RationalFunction], RationalFunction)
z) <- IndexedPolynomial
-> [RationalFunction]
-> IndexedPolynomial
-> Maybe ([RationalFunction], RationalFunction)
reduce IndexedPolynomial
x [] IndexedPolynomial
common = ([RationalFunction], RationalFunction)
z
  | Bool
otherwise = ([], RationalFunction
h) -- Should never happen, but a fallback if it does.
  where
    common :: IndexedPolynomial
common = IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
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 IndexedPolynomial
y (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate IndexedPolynomial
y
    (IndexedPolynomial
divisor, IndexedPolynomial
_) = IndexedPolynomial
y IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
common
    reduce :: IndexedPolynomial
-> [RationalFunction]
-> IndexedPolynomial
-> Maybe ([RationalFunction], RationalFunction)
reduce IndexedPolynomial
a [RationalFunction]
g IndexedPolynomial
d
      | IndexedPolynomial -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree IndexedPolynomial
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = do
          let d' :: IndexedPolynomial
d' = IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
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 IndexedPolynomial
d (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate IndexedPolynomial
d
          let (IndexedPolynomial
d'', IndexedPolynomial
_) = IndexedPolynomial
d IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
d'
          let (IndexedPolynomial
d''', IndexedPolynomial
_) = (IndexedPolynomial
divisor IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate IndexedPolynomial
d) IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
d
          (b, c) <- IndexedPolynomial
-> IndexedPolynomial
-> IndexedPolynomial
-> Maybe (IndexedPolynomial, IndexedPolynomial)
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 (-IndexedPolynomial
d''') IndexedPolynomial
d'' IndexedPolynomial
a
          let (b', _) = (differentiate b * divisor) `divide` d''
          let a' = IndexedPolynomial
c IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
- IndexedPolynomial
b'
          let g' = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
b IndexedPolynomial
d RationalFunction -> [RationalFunction] -> [RationalFunction]
forall a. a -> [a] -> [a]
: [RationalFunction]
g
          reduce a' g' d'
      | Bool
otherwise = ([RationalFunction], RationalFunction)
-> Maybe ([RationalFunction], RationalFunction)
forall a. a -> Maybe a
Just ([RationalFunction]
g, IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
a IndexedPolynomial
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 = toRationalFunction 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 :: RationalFunction
-> Maybe [(IndexedPolynomial, P Int IndexedPolynomial)]
rationalIntegralLogTerms (RationalFunction IndexedPolynomial
a IndexedPolynomial
d) = do
  -- For A/D, get the resultant and subresultant polynomial remainder sequence
  -- for D and (A - t * D').
  let sa :: P Int RationalFunction
sa = (Rational -> RationalFunction)
-> IndexedPolynomial -> P Int RationalFunction
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 -> RationalFunction
forall a. Fractional a => Rational -> a
fromRational IndexedPolynomial
a
  let sd :: P Int RationalFunction
sd = (Rational -> RationalFunction)
-> IndexedPolynomial -> P Int RationalFunction
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 -> RationalFunction
forall a. Fractional a => Rational -> a
fromRational IndexedPolynomial
d
  let t :: RationalFunction
t = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
RationalFunction (Int -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomial
1
  let (RationalFunction
resultant, [P Int RationalFunction]
prs) = P Int RationalFunction
-> P Int RationalFunction
-> (RationalFunction, [P Int RationalFunction])
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 RationalFunction
sd (P Int RationalFunction
 -> (RationalFunction, [P Int RationalFunction]))
-> P Int RationalFunction
-> (RationalFunction, [P Int RationalFunction])
forall a b. (a -> b) -> a -> b
$ P Int RationalFunction
sa P Int RationalFunction
-> P Int RationalFunction -> P Int RationalFunction
forall a. Num a => a -> a -> a
- RationalFunction
-> P Int RationalFunction -> P Int RationalFunction
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale RationalFunction
t (P Int RationalFunction -> P Int RationalFunction
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate P Int RationalFunction
sd)

  -- Turn rational functions into polynomials if possible.
  -- When the preconditions are satisfied, these should all be polynomials.
  sd' <- P Int RationalFunction -> Maybe (P Int IndexedPolynomial)
toPolyCoefficients P Int RationalFunction
sd
  resultant' <- toPoly resultant
  prs' <- toMaybeList $ map toPolyCoefficients prs :: Maybe [IndexedPolynomialWith IndexedPolynomial]

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

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

    derive ::
      IndexedPolynomial ->
      IndexedPolynomialWith IndexedPolynomial ->
      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
    derive :: IndexedPolynomial
-> P Int IndexedPolynomial
-> (IndexedPolynomial, P Int IndexedPolynomial)
derive IndexedPolynomial
q P Int IndexedPolynomial
s = (IndexedPolynomial
q, P Int IndexedPolynomial
s')
      where
        as :: [IndexedPolynomial]
as = IndexedPolynomial -> [IndexedPolynomial]
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 (IndexedPolynomial -> [IndexedPolynomial])
-> IndexedPolynomial -> [IndexedPolynomial]
forall a b. (a -> b) -> a -> b
$ P Int IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient P Int IndexedPolynomial
s
        s' :: P Int IndexedPolynomial
s' = (P Int IndexedPolynomial
 -> (Int, IndexedPolynomial) -> P Int IndexedPolynomial)
-> P Int IndexedPolynomial
-> [(Int, IndexedPolynomial)]
-> P Int IndexedPolynomial
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl P Int IndexedPolynomial
-> (Int, IndexedPolynomial) -> P Int IndexedPolynomial
forall {p :: * -> * -> *} {e} {p :: * -> * -> *} {b}.
(Polynomial p e IndexedPolynomial,
 Polynomial p e IndexedPolynomial, Integral b,
 Num (p e IndexedPolynomial)) =>
p e IndexedPolynomial
-> (b, IndexedPolynomial) -> p e IndexedPolynomial
scalePoly P Int IndexedPolynomial
s ([Int] -> [IndexedPolynomial] -> [(Int, IndexedPolynomial)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int
1 ..] :: [Int]) [IndexedPolynomial]
as)
          where
            scalePoly :: p e IndexedPolynomial
-> (b, IndexedPolynomial) -> p e IndexedPolynomial
scalePoly p e IndexedPolynomial
x (b
j, IndexedPolynomial
u) =
              Sum (p e IndexedPolynomial) -> p e IndexedPolynomial
forall a. Sum a -> a
getSum (Sum (p e IndexedPolynomial) -> p e IndexedPolynomial)
-> Sum (p e IndexedPolynomial) -> p e IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ (e -> IndexedPolynomial -> Sum (p e IndexedPolynomial))
-> p e IndexedPolynomial -> Sum (p e IndexedPolynomial)
forall m.
Monoid m =>
(e -> IndexedPolynomial -> m) -> p e IndexedPolynomial -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (IndexedPolynomial
-> e -> IndexedPolynomial -> Sum (p e IndexedPolynomial)
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 (IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (IndexedPolynomial -> IndexedPolynomial)
-> IndexedPolynomial -> IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
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 IndexedPolynomial
u IndexedPolynomial
q IndexedPolynomial -> b -> IndexedPolynomial
forall a b. (Num a, Integral b) => a -> b -> a
^ b
j)) p e IndexedPolynomial
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 -> IndexedPolynomial -> IndexedPolynomial -> Expression
complexLogTermToAtan Text
v IndexedPolynomial
a IndexedPolynomial
b
  | IndexedPolynomial
r IndexedPolynomial -> IndexedPolynomial -> Bool
forall a. Eq a => a -> a -> Bool
== IndexedPolynomial
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')
  | IndexedPolynomial -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree IndexedPolynomial
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< IndexedPolynomial -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree IndexedPolynomial
b = Text -> IndexedPolynomial -> IndexedPolynomial -> Expression
complexLogTermToAtan Text
v (-IndexedPolynomial
b) IndexedPolynomial
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 -> IndexedPolynomial -> IndexedPolynomial -> Expression
complexLogTermToAtan Text
v IndexedPolynomial
d IndexedPolynomial
c
  where
    (IndexedPolynomial
_, IndexedPolynomial
r) = IndexedPolynomial
a IndexedPolynomial
-> IndexedPolynomial -> (IndexedPolynomial, IndexedPolynomial)
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` IndexedPolynomial
b
    (IndexedPolynomial
d, IndexedPolynomial
c, IndexedPolynomial
g) = IndexedPolynomial
-> IndexedPolynomial
-> (IndexedPolynomial, IndexedPolynomial, IndexedPolynomial)
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 IndexedPolynomial
b (-IndexedPolynomial
a)
    a' :: Expression
a' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 IndexedPolynomial
a
    b' :: Expression
b' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 IndexedPolynomial
b
    g' :: Expression
g' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 IndexedPolynomial
g
    s' :: Expression
s' = Text -> (Rational -> Expression) -> IndexedPolynomial -> 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 (IndexedPolynomial -> Expression)
-> IndexedPolynomial -> Expression
forall a b. (a -> b) -> a -> b
$ IndexedPolynomial
a IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
d IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
+ IndexedPolynomial
b IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
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 :: (IndexedPolynomial, P Int IndexedPolynomial)
-> ((P Int IndexedPolynomial, P Int IndexedPolynomial),
    (IndexedPolynomialWith (P Int IndexedPolynomial),
     IndexedPolynomialWith (P Int IndexedPolynomial)))
complexLogTermToRealTerm (IndexedPolynomial
q, P Int IndexedPolynomial
s) = ((P Int IndexedPolynomial
qp, P Int IndexedPolynomial
qq), (IndexedPolynomialWith (P Int IndexedPolynomial)
sp, IndexedPolynomialWith (P Int IndexedPolynomial)
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 (P Int IndexedPolynomial)
q' = Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a. Sum a -> a
getSum (Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
 -> IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a b. (a -> b) -> a -> b
$ (Int
 -> P Int IndexedPolynomial
 -> Sum (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall m.
Monoid m =>
(Int -> P Int IndexedPolynomial -> m)
-> IndexedPolynomialWith (P Int IndexedPolynomial) -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> P Int IndexedPolynomial
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith (P Int IndexedPolynomial)
 -> Sum (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a. Sum a -> a
getSum (Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
 -> IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomial
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall m.
Monoid m =>
(Int -> Rational -> m) -> IndexedPolynomial -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> Rational
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
fromTerm IndexedPolynomial
q
      where
        fromTerm :: Int -> Rational -> Sum (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
        fromTerm :: Int
-> Rational
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
fromTerm Int
e Rational
c = IndexedPolynomialWith (P Int IndexedPolynomial)
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. a -> Sum a
Sum (IndexedPolynomialWith (P Int IndexedPolynomial)
 -> Sum (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> Sum (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (P Int IndexedPolynomial)
c' IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith (P Int IndexedPolynomial)
u IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith (P Int IndexedPolynomial)
i IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a. Num a => a -> a -> a
* IndexedPolynomialWith (P Int IndexedPolynomial)
v) IndexedPolynomialWith (P Int IndexedPolynomial)
-> Int -> IndexedPolynomialWith (P Int IndexedPolynomial)
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith (P Int IndexedPolynomial)
c' = P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> P Int IndexedPolynomial -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational -> IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c IndexedPolynomial
1) P Int IndexedPolynomial
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1
        i :: IndexedPolynomialWith (P Int IndexedPolynomial)
i = Int -> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        u :: IndexedPolynomialWith (P Int IndexedPolynomial)
u = P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1
        v :: IndexedPolynomialWith (P Int IndexedPolynomial)
v = P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> P Int IndexedPolynomial -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) P Int IndexedPolynomial
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1
    -- q' == qp + i * qq
    (P Int IndexedPolynomial
qp, P Int IndexedPolynomial
qq) = (IndexedPolynomialWith (P Int IndexedPolynomial)
-> Int -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (P Int IndexedPolynomial)
q' Int
0, IndexedPolynomialWith (P Int IndexedPolynomial)
-> Int -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (P Int IndexedPolynomial)
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 (P Int IndexedPolynomial))
s' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomialWith (P Int IndexedPolynomial)
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall m.
Monoid m =>
(Int -> IndexedPolynomialWith (P Int IndexedPolynomial) -> m)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith
   (IndexedPolynomialWith (P Int IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomial
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> P Int IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall m.
Monoid m =>
(Int -> IndexedPolynomial -> m) -> P Int IndexedPolynomial -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
fromTerm P Int IndexedPolynomial
s
      where
        fromTerm :: Int -> IndexedPolynomial -> Sum (IndexedPolynomialWith (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
        fromTerm :: Int
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
fromTerm Int
e IndexedPolynomial
c = IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (P Int IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
c' IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
x IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
c' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall m.
Monoid m =>
(Int -> Rational -> m) -> IndexedPolynomial -> 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 IndexedPolynomial)))
forall {b}.
Integral b =>
b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
fromCoefficient IndexedPolynomial
c
            fromCoefficient :: b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
fromCoefficient b
e' Rational
c'' = IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (P Int IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (P Int IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (P Int IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
c''' IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
u IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
i IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
v) IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> b
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall a b. (Num a, Integral b) => a -> b -> a
^ b
e'
              where
                c''' :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
c''' = IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> P Int IndexedPolynomial -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational -> IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c'' IndexedPolynomial
1) P Int IndexedPolynomial
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
1
        i :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
i = Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        x :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
x = IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
1
        u :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
u = IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
1
        v :: IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
v = IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (P Int IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (P Int IndexedPolynomial
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> P Int IndexedPolynomial -> P Int IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) P Int IndexedPolynomial
1) IndexedPolynomialWith (P Int IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
1
    -- s' = sp + i * sq
    (IndexedPolynomialWith (P Int IndexedPolynomial)
sp, IndexedPolynomialWith (P Int IndexedPolynomial)
sq) = (IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Int -> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
s' Int
0, IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
-> Int -> IndexedPolynomialWith (P Int IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (P Int IndexedPolynomial))
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
-> (IndexedPolynomial, P Int IndexedPolynomial) -> Maybe Expression
complexLogTermToRealExpression Text
v (IndexedPolynomial
r, P Int IndexedPolynomial
s)
  | (Just [(Rational, Rational)]
xys) <- P Int IndexedPolynomial
-> P Int IndexedPolynomial -> Maybe [(Rational, Rational)]
solveBivariatePolynomials P Int IndexedPolynomial
p P Int IndexedPolynomial
q,
    (Just [Expression]
h) <- [(Rational, Rational)] -> Maybe [Expression]
f [(Rational, Rational)]
xys,
    (Just [Rational]
zs) <- Maybe [Expression] -> Maybe [Rational]
toRationalList (IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
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
    ((P Int IndexedPolynomial
p, P Int IndexedPolynomial
q), (IndexedPolynomialWith (P Int IndexedPolynomial)
a, IndexedPolynomialWith (P Int IndexedPolynomial)
b)) = (IndexedPolynomial, P Int IndexedPolynomial)
-> ((P Int IndexedPolynomial, P Int IndexedPolynomial),
    (IndexedPolynomialWith (P Int IndexedPolynomial),
     IndexedPolynomialWith (P Int IndexedPolynomial)))
complexLogTermToRealTerm (IndexedPolynomial
r, P Int IndexedPolynomial
s)

    f :: [(Rational, Rational)] -> Maybe [Expression]
    f :: [(Rational, Rational)] -> Maybe [Expression]
f [(Rational, Rational)]
xys = [Maybe Expression] -> Maybe [Expression]
forall a. [Maybe a] -> Maybe [a]
toMaybeList ([Maybe Expression] -> Maybe [Expression])
-> [Maybe Expression] -> Maybe [Expression]
forall a b. (a -> b) -> a -> b
$ do
      (x, 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'' = (IndexedPolynomial -> Expression)
-> P Int IndexedPolynomial -> 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) -> IndexedPolynomial -> 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' = (P Int IndexedPolynomial -> Expression)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> 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)
-> (P Int IndexedPolynomial -> P Int Expression)
-> P Int IndexedPolynomial
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. P Int IndexedPolynomial -> P Int Expression
flatten'') -- u-polynomials into Expressions.
      let 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 (P Int IndexedPolynomial)
    -> P Int Expression)
-> IndexedPolynomialWith (P Int IndexedPolynomial)
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith (P Int IndexedPolynomial) -> P Int Expression
flatten' -- x-polynomials into Expressions.
      -- a and b flattened into Expressions.
      let a' = IndexedPolynomialWith (P Int IndexedPolynomial) -> Expression
flatten IndexedPolynomialWith (P Int IndexedPolynomial)
a
      let b' = IndexedPolynomialWith (P Int IndexedPolynomial) -> Expression
flatten IndexedPolynomialWith (P Int IndexedPolynomial)
b
      -- a and b flattened into x-polynomials with rational number coefficients.
      return $ do
        a'' <- convertCoefficients $ flatten' a
        b'' <- convertCoefficients $ flatten' b
        return $ fromRational x * log (a' * a' + b' * b') + fromRational y * complexLogTermToAtan v a'' 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
      z <- t Rational
zs
      let s' = (IndexedPolynomial -> Expression)
-> P Int IndexedPolynomial -> 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) -> IndexedPolynomial -> 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) P Int IndexedPolynomial
s
      return $ fromRational z * Log' (toExpression v toSymbolicCoefficient 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 IndexedPolynomial
convertCoefficients P Int Expression
x = [IndexedPolynomial] -> IndexedPolynomial
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([IndexedPolynomial] -> IndexedPolynomial)
-> ([(Int, Rational)] -> [IndexedPolynomial])
-> [(Int, Rational)]
-> IndexedPolynomial
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Rational) -> IndexedPolynomial)
-> [(Int, Rational)] -> [IndexedPolynomial]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
e, Rational
c) -> Rational -> IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale Rational
c (Int -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
e)) ([(Int, Rational)] -> IndexedPolynomial)
-> Maybe [(Int, Rational)] -> Maybe IndexedPolynomial
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Maybe (Int, Rational)] -> Maybe [(Int, Rational)]
forall a. [Maybe a] -> Maybe [a]
toMaybeList ((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

-- | 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 :: P Int IndexedPolynomial
-> P Int IndexedPolynomial -> Maybe [(Rational, Rational)]
solveBivariatePolynomials P Int IndexedPolynomial
p P Int IndexedPolynomial
q = do
  let p' :: P Int RationalFunction
p' = P Int IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients P Int IndexedPolynomial
p
  let q' :: P Int RationalFunction
q' = P Int IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients P Int IndexedPolynomial
q
  resultant <- RationalFunction -> Maybe IndexedPolynomial
toPoly (RationalFunction -> Maybe IndexedPolynomial)
-> RationalFunction -> Maybe IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ (RationalFunction, [P Int RationalFunction]) -> RationalFunction
forall a b. (a, b) -> a
fst ((RationalFunction, [P Int RationalFunction]) -> RationalFunction)
-> (RationalFunction, [P Int RationalFunction]) -> RationalFunction
forall a b. (a -> b) -> a -> b
$ P Int RationalFunction
-> P Int RationalFunction
-> (RationalFunction, [P Int RationalFunction])
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 RationalFunction
p' P Int RationalFunction
q'
  vs' <- solve resultant
  vs <- toMaybeList $ map (convert . simplify) vs'
  concat <$> toMaybeList (map solveForU vs)
  where
    toRationalFunctionCoefficients :: P Int IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients = (IndexedPolynomial -> RationalFunction)
-> P Int IndexedPolynomial -> P Int RationalFunction
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 (IndexedPolynomial -> IndexedPolynomial -> RationalFunction
`toRationalFunction` IndexedPolynomial
1)

    -- 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
      | IndexedPolynomial
0 <- IndexedPolynomial
p' = do
          -- Any u will make p'=0 true, so we only need to solve p'.
          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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
q'
          map (,v) <$> toMaybeList u
      | IndexedPolynomial
0 <- IndexedPolynomial
q' = do
          -- Any u will make q'=0 true, so we only need to solve p'.
          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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
p'
          map (,v) <$> toMaybeList u
      | Bool
otherwise = do
          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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
p'
          uq <- map (convert . simplify) <$> solve q'
          up' <- toMaybeList up
          uq' <- toMaybeList uq
          return $ map (,v) $ up' `intersect` uq'
      where
        p' :: IndexedPolynomial
p' = (IndexedPolynomial -> Rational)
-> P Int IndexedPolynomial -> IndexedPolynomial
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)
-> (IndexedPolynomial -> Sum Rational)
-> IndexedPolynomial
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Rational -> Sum Rational)
-> IndexedPolynomial -> Sum Rational
forall m.
Monoid m =>
(Int -> Rational -> m) -> IndexedPolynomial -> 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)) P Int IndexedPolynomial
p
        q' :: IndexedPolynomial
q' = (IndexedPolynomial -> Rational)
-> P Int IndexedPolynomial -> IndexedPolynomial
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)
-> (IndexedPolynomial -> Sum Rational)
-> IndexedPolynomial
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Rational -> Sum Rational)
-> IndexedPolynomial -> Sum Rational
forall m.
Monoid m =>
(Int -> Rational -> m) -> IndexedPolynomial -> 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)) P Int IndexedPolynomial
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

-- | Turn the rational function into a polynomial if possible.
toPoly :: RationalFunction -> Maybe IndexedPolynomial
toPoly :: RationalFunction -> Maybe IndexedPolynomial
toPoly (RationalFunction IndexedPolynomial
p IndexedPolynomial
q)
  | IndexedPolynomial -> Int
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree IndexedPolynomial
q Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0, IndexedPolynomial
q IndexedPolynomial -> IndexedPolynomial -> Bool
forall a. Eq a => a -> a -> Bool
/= IndexedPolynomial
0 = IndexedPolynomial -> Maybe IndexedPolynomial
forall a. a -> Maybe a
Just IndexedPolynomial
p'
  | Bool
otherwise = Maybe IndexedPolynomial
forall a. Maybe a
Nothing
  where
    p' :: IndexedPolynomial
p' = Rational -> IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ IndexedPolynomial -> Rational
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient IndexedPolynomial
q) IndexedPolynomial
p

-- | Turn the rational function coefficients into polynomial coefficients if possible.
toPolyCoefficients ::
  IndexedPolynomialWith RationalFunction ->
  Maybe (IndexedPolynomialWith IndexedPolynomial)
toPolyCoefficients :: P Int RationalFunction -> Maybe (P Int IndexedPolynomial)
toPolyCoefficients P Int RationalFunction
p = [(Int, Maybe IndexedPolynomial)] -> Maybe (P Int IndexedPolynomial)
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Num (p e c)) =>
[(e, Maybe c)] -> Maybe (p e c)
reconstruct [(Int, Maybe IndexedPolynomial)]
terms
  where
    terms :: [(Int, Maybe IndexedPolynomial)]
terms = (Int -> RationalFunction -> [(Int, Maybe IndexedPolynomial)])
-> P Int RationalFunction -> [(Int, Maybe IndexedPolynomial)]
forall m.
Monoid m =>
(Int -> RationalFunction -> m) -> P Int RationalFunction -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms (\Int
e RationalFunction
c -> [(Int
e, RationalFunction -> Maybe IndexedPolynomial
toPoly RationalFunction
c)]) P Int RationalFunction
p
    reconstruct :: [(e, Maybe c)] -> Maybe (p e c)
reconstruct [] = p e c -> Maybe (p e c)
forall a. a -> Maybe a
Just p e c
0
    reconstruct ((e
_, Maybe c
Nothing) : [(e, Maybe c)]
_) = Maybe (p e c)
forall a. Maybe a
Nothing
    reconstruct ((e
e, Just c
c) : [(e, Maybe c)]
xs)
      | (Just p e c
p') <- [(e, Maybe c)] -> Maybe (p e c)
reconstruct [(e, Maybe c)]
xs = p e c -> Maybe (p e c)
forall a. a -> Maybe a
Just (p e c -> Maybe (p e c)) -> p e c -> Maybe (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c (e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
e) p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
p'
      | Bool
otherwise = Maybe (p e c)
forall a. Maybe a
Nothing

-- | If there are any nothings, then turn the list into nothing.
-- Otherwise, turn it into the list of just the elements.
toMaybeList :: [Maybe a] -> Maybe [a]
toMaybeList :: forall a. [Maybe a] -> Maybe [a]
toMaybeList [] = [a] -> Maybe [a]
forall a. a -> Maybe a
Just []
toMaybeList (Maybe a
Nothing : [Maybe a]
_) = Maybe [a]
forall a. Maybe a
Nothing
toMaybeList (Just a
x : [Maybe a]
xs)
  | (Just [a]
xs') <- [Maybe a] -> Maybe [a]
forall a. [Maybe a] -> Maybe [a]
toMaybeList [Maybe a]
xs = [a] -> Maybe [a]
forall a. a -> Maybe a
Just (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xs')
  | Bool
otherwise = Maybe [a]
forall a. Maybe a
Nothing