-- |
-- 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.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.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)\]
--
-- 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 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, IndexedPolynomialWith IndexedPolynomial)]
logTerms = RationalFunction
-> Maybe
     [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
rationalIntegralLogTerms RationalFunction
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 [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
terms) <- Maybe
  [(IndexedPolynomial, IndexedPolynomialWith 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, IndexedPolynomialWith IndexedPolynomial)
 -> Maybe Expression)
-> [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
-> [Maybe Expression]
forall a b. (a -> b) -> [a] -> [b]
map (Text
-> (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
-> Maybe Expression
complexLogTermToRealExpression Text
v) [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
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 [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
terms) <- Maybe
  [(IndexedPolynomial, IndexedPolynomialWith 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, IndexedPolynomialWith IndexedPolynomial)
 -> Maybe Expression)
-> [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]
-> [Maybe Expression]
forall a b. (a -> b) -> [a] -> [b]
map (Text
-> (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
-> Maybe Expression
complexLogTermToComplexExpression Text
v) [(IndexedPolynomial, IndexedPolynomialWith 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
          (IndexedPolynomial
b, IndexedPolynomial
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 (IndexedPolynomial
b', IndexedPolynomial
_) = (IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate IndexedPolynomial
b IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
* IndexedPolynomial
divisor) 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 a' :: IndexedPolynomial
a' = IndexedPolynomial
c IndexedPolynomial -> IndexedPolynomial -> IndexedPolynomial
forall a. Num a => a -> a -> a
- IndexedPolynomial
b'
          let g' :: [RationalFunction]
g' = IndexedPolynomial -> IndexedPolynomial -> RationalFunction
toRationalFunction IndexedPolynomial
b IndexedPolynomial
d RationalFunction -> [RationalFunction] -> [RationalFunction]
forall a. a -> [a] -> [a]
: [RationalFunction]
g
          IndexedPolynomial
-> [RationalFunction]
-> IndexedPolynomial
-> Maybe ([RationalFunction], RationalFunction)
reduce IndexedPolynomial
a' [RationalFunction]
g' IndexedPolynomial
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, IndexedPolynomialWith 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.
  IndexedPolynomialWith IndexedPolynomial
sd' <- (RationalFunction -> Maybe IndexedPolynomial)
-> P Int RationalFunction
-> Maybe (IndexedPolynomialWith IndexedPolynomial)
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 RationalFunction -> Maybe IndexedPolynomial
toPoly P Int RationalFunction
sd
  IndexedPolynomial
resultant' <- RationalFunction -> Maybe IndexedPolynomial
toPoly RationalFunction
resultant
  [IndexedPolynomialWith IndexedPolynomial]
prs' <- [Maybe (IndexedPolynomialWith IndexedPolynomial)]
-> Maybe [IndexedPolynomialWith IndexedPolynomial]
forall a. [Maybe a] -> Maybe [a]
toMaybeList ([Maybe (IndexedPolynomialWith IndexedPolynomial)]
 -> Maybe [IndexedPolynomialWith IndexedPolynomial])
-> [Maybe (IndexedPolynomialWith IndexedPolynomial)]
-> Maybe [IndexedPolynomialWith IndexedPolynomial]
forall a b. (a -> b) -> a -> b
$ (P Int RationalFunction
 -> Maybe (IndexedPolynomialWith IndexedPolynomial))
-> [P Int RationalFunction]
-> [Maybe (IndexedPolynomialWith IndexedPolynomial)]
forall a b. (a -> b) -> [a] -> [b]
map ((RationalFunction -> Maybe IndexedPolynomial)
-> P Int RationalFunction
-> Maybe (IndexedPolynomialWith IndexedPolynomial)
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 RationalFunction -> Maybe IndexedPolynomial
toPoly) [P Int RationalFunction]
prs :: Maybe [IndexedPolynomialWith IndexedPolynomial]

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

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

    derive ::
      IndexedPolynomial ->
      IndexedPolynomialWith IndexedPolynomial ->
      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
    derive :: IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
derive IndexedPolynomial
q IndexedPolynomialWith IndexedPolynomial
s = (IndexedPolynomial
q, IndexedPolynomialWith 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
$ IndexedPolynomialWith IndexedPolynomial -> IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient IndexedPolynomialWith IndexedPolynomial
s
        s' :: IndexedPolynomialWith IndexedPolynomial
s' = (IndexedPolynomialWith IndexedPolynomial
 -> (Int, IndexedPolynomial)
 -> IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith IndexedPolynomial
-> [(Int, IndexedPolynomial)]
-> IndexedPolynomialWith IndexedPolynomial
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl IndexedPolynomialWith IndexedPolynomial
-> (Int, IndexedPolynomial)
-> IndexedPolynomialWith 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 IndexedPolynomialWith 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, IndexedPolynomialWith IndexedPolynomial)
-> ((IndexedPolynomialWith IndexedPolynomial,
     IndexedPolynomialWith IndexedPolynomial),
    (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial),
     IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
complexLogTermToRealTerm (IndexedPolynomial
q, IndexedPolynomialWith IndexedPolynomial
s) = ((IndexedPolynomialWith IndexedPolynomial
qp, IndexedPolynomialWith IndexedPolynomial
qq), (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
sp, IndexedPolynomialWith (IndexedPolynomialWith 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 (IndexedPolynomialWith IndexedPolynomial)
q' = Sum
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
 -> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomialWith IndexedPolynomial
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall m.
Monoid m =>
(Int -> IndexedPolynomialWith IndexedPolynomial -> m)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomialWith IndexedPolynomial
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ Sum
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
 -> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith 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 IndexedPolynomial))
fromTerm IndexedPolynomial
q
      where
        fromTerm :: Int -> Rational -> Sum (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
        fromTerm :: Int
-> Rational
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
fromTerm Int
e Rational
c = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. a -> Sum a
Sum (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
 -> Sum
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
c' IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
u IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
i IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a. Num a => a -> a -> a
* IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
v) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
c' = IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith 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) IndexedPolynomialWith IndexedPolynomial
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1
        i :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
i = Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        u :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
u = IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomialWith IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1
        v :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
v = IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith 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) IndexedPolynomialWith IndexedPolynomial
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1
    -- q' == qp + i * qq
    (IndexedPolynomialWith IndexedPolynomial
qp, IndexedPolynomialWith IndexedPolynomial
qq) = (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Int -> IndexedPolynomialWith IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
q' Int
0, IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Int -> IndexedPolynomialWith IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith (IndexedPolynomialWith 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 (IndexedPolynomialWith IndexedPolynomial))
s' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall m.
Monoid m =>
(Int
 -> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
 -> m)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a.
(Eq a, Num a) =>
Int -> a -> Sum (IndexedPolynomialWith a)
reduceImaginary (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> IndexedPolynomial
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomialWith IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall m.
Monoid m =>
(Int -> IndexedPolynomial -> m)
-> IndexedPolynomialWith 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 (IndexedPolynomialWith IndexedPolynomial)))
fromTerm IndexedPolynomialWith IndexedPolynomial
s
      where
        fromTerm :: Int -> IndexedPolynomial -> Sum (IndexedPolynomialWith (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
        fromTerm :: Int
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
fromTerm Int
e IndexedPolynomial
c = IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
c' IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
x IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e
          where
            c' :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
c' = Sum
  (IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Sum a -> a
getSum (Sum
   (IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
 -> IndexedPolynomialWith
      (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (a -> b) -> a -> b
$ (Int
 -> Rational
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomial
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith 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 (IndexedPolynomialWith IndexedPolynomial)))
forall {b}.
Integral b =>
b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
fromCoefficient IndexedPolynomial
c
            fromCoefficient :: b
-> Rational
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
fromCoefficient b
e' Rational
c'' = IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a. a -> Sum a
Sum (IndexedPolynomialWith
   (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
 -> Sum
      (IndexedPolynomialWith
         (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Sum
     (IndexedPolynomialWith
        (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
c''' IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Num a => a -> a -> a
* (IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
u IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Num a => a -> a -> a
+ IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
i IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a. Num a => a -> a -> a
* IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
v) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> b
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall a b. (Num a, Integral b) => a -> b -> a
^ b
e'
              where
                c''' :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
c''' = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith 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) IndexedPolynomialWith IndexedPolynomial
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
1
        i :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
i = Int
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1
        x :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
x = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
1
        u :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
u = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Int -> IndexedPolynomialWith IndexedPolynomial
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power Int
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
1
        v :: IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
v = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> IndexedPolynomialWith
     (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith 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) IndexedPolynomialWith IndexedPolynomial
1) IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
1) IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
1
    -- s' = sp + i * sq
    (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
sp, IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
sq) = (IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
s' Int
0, IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))
-> Int
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e -> c
coefficient IndexedPolynomialWith
  (IndexedPolynomialWith (IndexedPolynomialWith 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, IndexedPolynomialWith IndexedPolynomial)
-> Maybe Expression
complexLogTermToRealExpression Text
v (IndexedPolynomial
r, IndexedPolynomialWith IndexedPolynomial
s)
  | (Just [(Rational, Rational)]
xys) <- IndexedPolynomialWith IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> Maybe [(Rational, Rational)]
solveBivariatePolynomials IndexedPolynomialWith IndexedPolynomial
p IndexedPolynomialWith 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
    ((IndexedPolynomialWith IndexedPolynomial
p, IndexedPolynomialWith IndexedPolynomial
q), (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
a, IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
b)) = (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
-> ((IndexedPolynomialWith IndexedPolynomial,
     IndexedPolynomialWith IndexedPolynomial),
    (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial),
     IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))
complexLogTermToRealTerm (IndexedPolynomial
r, IndexedPolynomialWith 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
      (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 IndexedPolynomial -> P Int Expression
flatten'' = (IndexedPolynomial -> Expression)
-> IndexedPolynomialWith 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' :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> P Int Expression
flatten' = (IndexedPolynomialWith IndexedPolynomial -> Expression)
-> IndexedPolynomialWith (IndexedPolynomialWith 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)
-> (IndexedPolynomialWith IndexedPolynomial -> P Int Expression)
-> IndexedPolynomialWith IndexedPolynomial
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith IndexedPolynomial -> P Int Expression
flatten'') -- u-polynomials into Expressions.
      let flatten :: IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> 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 IndexedPolynomial)
    -> P Int Expression)
-> IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Expression
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> P Int Expression
flatten' -- x-polynomials into Expressions.
      -- a and b flattened into Expressions.
      let a' :: Expression
a' = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Expression
flatten IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
a
      let b' :: Expression
b' = IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> Expression
flatten IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
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
        IndexedPolynomial
a'' <- P Int Expression -> Maybe IndexedPolynomial
convertCoefficients (P Int Expression -> Maybe IndexedPolynomial)
-> P Int Expression -> Maybe IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> P Int Expression
flatten' IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
a
        IndexedPolynomial
b'' <- P Int Expression -> Maybe IndexedPolynomial
convertCoefficients (P Int Expression -> Maybe IndexedPolynomial)
-> P Int Expression -> Maybe IndexedPolynomial
forall a b. (a -> b) -> a -> b
$ IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
-> P Int Expression
flatten' IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)
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 -> IndexedPolynomial -> IndexedPolynomial -> Expression
complexLogTermToAtan Text
v IndexedPolynomial
a'' IndexedPolynomial
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' = (IndexedPolynomial -> Expression)
-> IndexedPolynomialWith 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) IndexedPolynomialWith IndexedPolynomial
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 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

-- | 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
-> (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)
-> Maybe Expression
complexLogTermToComplexExpression Text
v (IndexedPolynomial
q, IndexedPolynomialWith IndexedPolynomial
s) = do
  [Expression]
as <- IndexedPolynomial -> Maybe [Expression]
complexSolve IndexedPolynomial
q
  let terms :: [Expression]
terms = do
        Expression
a <- [Expression]
as
        let s' :: P Int Expression
s' = (IndexedPolynomial -> Expression)
-> IndexedPolynomialWith 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 -> IndexedPolynomial -> Expression
forall {p :: * -> * -> *} {a} {a}.
(Polynomial p a Rational, Floating a) =>
a -> p a Rational -> a
collapse Expression
a) IndexedPolynomialWith IndexedPolynomial
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 IndexedPolynomial
-> IndexedPolynomialWith IndexedPolynomial
-> Maybe [(Rational, Rational)]
solveBivariatePolynomials IndexedPolynomialWith IndexedPolynomial
p IndexedPolynomialWith IndexedPolynomial
q = do
  let p' :: P Int RationalFunction
p' = IndexedPolynomialWith IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients IndexedPolynomialWith IndexedPolynomial
p
  let q' :: P Int RationalFunction
q' = IndexedPolynomialWith IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients IndexedPolynomialWith IndexedPolynomial
q
  IndexedPolynomial
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'
  [Expression]
vs' <- IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
resultant
  [Rational]
vs <- [Maybe Rational] -> Maybe [Rational]
forall a. [Maybe a] -> Maybe [a]
toMaybeList ([Maybe Rational] -> Maybe [Rational])
-> [Maybe Rational] -> Maybe [Rational]
forall a b. (a -> b) -> a -> b
$ (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]
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
<$> [Maybe [(Rational, Rational)]] -> Maybe [[(Rational, Rational)]]
forall a. [Maybe a] -> Maybe [a]
toMaybeList ((Rational -> Maybe [(Rational, Rational)])
-> [Rational] -> [Maybe [(Rational, Rational)]]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Maybe [(Rational, Rational)]
solveForU [Rational]
vs)
  where
    toRationalFunctionCoefficients :: IndexedPolynomialWith IndexedPolynomial -> P Int RationalFunction
toRationalFunctionCoefficients = (IndexedPolynomial -> RationalFunction)
-> IndexedPolynomialWith 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'.
          [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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
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 a. [Maybe a] -> Maybe [a]
toMaybeList [Maybe Rational]
u
      | IndexedPolynomial
0 <- IndexedPolynomial
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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
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 a. [Maybe a] -> Maybe [a]
toMaybeList [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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
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
<$> IndexedPolynomial -> Maybe [Expression]
solve IndexedPolynomial
q'
          [Rational]
up' <- [Maybe Rational] -> Maybe [Rational]
forall a. [Maybe a] -> Maybe [a]
toMaybeList [Maybe Rational]
up
          [Rational]
uq' <- [Maybe Rational] -> Maybe [Rational]
forall a. [Maybe a] -> Maybe [a]
toMaybeList [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' :: IndexedPolynomial
p' = (IndexedPolynomial -> Rational)
-> IndexedPolynomialWith 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)) IndexedPolynomialWith IndexedPolynomial
p
        q' :: IndexedPolynomial
q' = (IndexedPolynomial -> Rational)
-> IndexedPolynomialWith 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)) IndexedPolynomialWith 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

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