{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE PatternSynonyms #-}

-- |
-- Module: Symtegration.Polynomial.Rational
-- Description: Representation for rational functions.
-- Copyright: Copyright 2025 Yoo Chung
-- License: Apache-2.0
-- Maintainer: dev@chungyc.org
--
-- Provides a representation for rational functions,
-- which are functions of ratios for two polynomials.
-- Note that these are /not/ functions of rational numbers.
module Symtegration.Polynomial.Rational
  ( Function,
    fromPolynomial,
    fromPolynomials,
    toPolynomial,
    pattern Function,
  )
where

import Control.DeepSeq (NFData)
import Data.Text (unpack)
import GHC.Generics (Generic)
import Symtegration.Polynomial
import TextShow

-- $setup
-- >>> import Symtegration.Polynomial
-- >>> import Symtegration.Polynomial.Indexed

-- | Represents a rational function.
--
-- Values can be constructed with 'fromPolynomial' and 'fromPolynomials',
-- while the numerator and denominator can be obtained by pattern matching with 'Function'.
-- Internally, common factors will be canceled out and the denominator will be monic.
--
-- >>> fromPolynomial $ power 1 :: Function IndexedPolynomial
-- Function (x) (1)
-- >>> fromPolynomials (power 2 + 1) (power 3 + 1) :: Function IndexedPolynomial
-- Function (x^2 + 1) (x^3 + 1)
-- >>> let (Function p q) = fromPolynomials (power 2) (power 1 + 1) :: Function IndexedPolynomial
-- >>> (p, q)
-- (x^2,x + 1)
-- >>> fromPolynomials (4 * power 2 - 4) (2 * power 1 - 2 :: IndexedPolynomial)
-- Function (2x + 2) (1)
--
-- The type is an instance of the 'Fractional' type class, so values can be added,
-- subtracted, multiplied, and divided.  Their values can also be given as
-- integer literals as well.
--
-- >>> let p = fromPolynomials (power 2 + 1) (power 3 + 1) :: Function IndexedPolynomial
-- >>> let q = fromPolynomials (power 1) (power 2 - 1) :: Function IndexedPolynomial
-- >>> p + q
-- Function (2x^3 + (-2)x^2 + 2x + (-1)) (x^4 + (-1)x^3 + x + (-1))
-- >>> p * q
-- Function (x^3 + x) (x^5 + (-1)x^3 + x^2 + (-1))
-- >>> 0 :: Function IndexedPolynomial
-- Function (0) (1)
-- >>> 5 :: Function IndexedPolynomial
-- Function (5) (1)
--
-- One could import this module with an alias so that @Rational.Function@ can be
-- used as the name of the type for rational functions.
--
-- >>> import Symtegration.Polynomial.Rational as Rational
-- >>> fromPolynomials (power 1) (power 2 + 1) :: Rational.Function IndexedPolynomial
-- Function (x) (x^2 + 1)
data Function a = F a a deriving ((forall x. Function a -> Rep (Function a) x)
-> (forall x. Rep (Function a) x -> Function a)
-> Generic (Function a)
forall x. Rep (Function a) x -> Function a
forall x. Function a -> Rep (Function a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (Function a) x -> Function a
forall a x. Function a -> Rep (Function a) x
$cfrom :: forall a x. Function a -> Rep (Function a) x
from :: forall x. Function a -> Rep (Function a) x
$cto :: forall a x. Rep (Function a) x -> Function a
to :: forall x. Rep (Function a) x -> Function a
Generic, Function a -> ()
(Function a -> ()) -> NFData (Function a)
forall a. NFData a => Function a -> ()
forall a. (a -> ()) -> NFData a
$crnf :: forall a. NFData a => Function a -> ()
rnf :: Function a -> ()
NFData)

{-# COMPLETE Function #-}

-- | Pattern synonym for matching the numerator and denominator of a rational function.
-- This would be equivalent to a data constructor,
-- except it cannot be used for constructing a rational function.
--
-- >>> let (Function p q) = fromPolynomials (power 2) (power 1 + 1) :: Function IndexedPolynomial
-- >>> p
-- x^2
-- >>> q
-- x + 1
pattern Function :: a -> a -> Function a
pattern $mFunction :: forall {r} {a}. Function a -> (a -> a -> r) -> ((# #) -> r) -> r
Function x y <- F x y

instance (Show a, TextShow a) => Show (Function a) where
  show :: Function a -> String
show = Text -> String
unpack (Text -> String) -> (Function a -> Text) -> Function a -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Function a -> Text
forall a. TextShow a => a -> Text
showt

instance (TextShow a) => TextShow (Function a) where
  showb :: Function a -> Builder
showb (F a
x a
y) = Builder
"Function (" Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> a -> Builder
forall a. TextShow a => a -> Builder
showb a
x Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Builder
") (" Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> a -> Builder
forall a. TextShow a => a -> Builder
showb a
y Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Builder
")"

instance (Eq a, Num a) => Eq (Function a) where
  (F a
x a
0) == :: Function a -> Function a -> Bool
== (F a
y a
0) = a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
y
  (F a
x a
y) == (F a
u a
v) = a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
v a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
y a -> a -> a
forall a. Num a => a -> a -> a
* a
u

instance (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) => Num (Function (p e c)) where
  (F p e c
x p e c
y) + :: Function (p e c) -> Function (p e c) -> Function (p e c)
+ (F p e c
u p e c
v) = p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials (p e c
x p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
u p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
y) (p e c
y p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v)
  (F p e c
x p e c
y) - :: Function (p e c) -> Function (p e c) -> Function (p e c)
- (F p e c
u p e c
v) = p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials (p e c
x p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c
u p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
y) (p e c
y p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v)
  (F p e c
x p e c
y) * :: Function (p e c) -> Function (p e c) -> Function (p e c)
* (F p e c
u p e c
v) = p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials (p e c
x p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
u) (p e c
y p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v)
  abs :: Function (p e c) -> Function (p e c)
abs = Function (p e c) -> Function (p e c)
forall a. a -> a
id
  signum :: Function (p e c) -> Function (p e c)
signum = Function (p e c) -> Function (p e c) -> Function (p e c)
forall a b. a -> b -> a
const Function (p e c)
1
  fromInteger :: Integer -> Function (p e c)
fromInteger Integer
n = p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c)) -> p e c -> Function (p e c)
forall a b. (a -> b) -> a -> b
$ Integer -> p e c
forall a. Num a => Integer -> a
fromInteger Integer
n

instance (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) => Fractional (Function (p e c)) where
  (F p e c
x p e c
y) / :: Function (p e c) -> Function (p e c) -> Function (p e c)
/ (F p e c
u p e c
v) = p e c -> p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials (p e c
x p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
v) (p e c
y p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
u)
  fromRational :: Rational -> Function (p e c)
fromRational Rational
x = p e c -> Function (p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial (p e c -> Function (p e c)) -> p e c -> Function (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (Rational -> c
forall a. Fractional a => Rational -> a
fromRational Rational
x) p e c
1

-- | Returns an equivalent rational function from a given polynomial.
--
-- >>> fromPolynomial (power 2 + 1 :: IndexedPolynomial)
-- Function (x^2 + 1) (1)
fromPolynomial :: (Polynomial p e c, Num (p e c)) => p e c -> Function (p e c)
fromPolynomial :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c)) =>
p e c -> Function (p e c)
fromPolynomial p e c
x = p e c -> p e c -> Function (p e c)
forall a. a -> a -> Function a
F p e c
x p e c
1

-- | Returns a rational function with the ratio of two polynomials.
--
-- >>> fromPolynomials (power 1) (power 2 + 1 :: IndexedPolynomial)
-- Function (x) (x^2 + 1)
fromPolynomials ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
  -- | Numerator \(a\).
  p e c ->
  -- | Denominator \(d\).
  p e c ->
  -- | Rational function \(\frac{a}{d}\).
  Function (p e c)
fromPolynomials :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> Function (p e c)
fromPolynomials p e c
x p e c
y = p e c -> p e c -> Function (p e c)
forall a. a -> a -> Function a
F p e c
x'' p e c
y''
  where
    -- Cancel out common factors.
    g :: p e c
g = p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
x p e c
y
    (p e c
x', p e c
_) = p e c
x 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
g
    (p e c
y', p e c
_) = p e c
y p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
g

    -- Make the denominator monic.
    (p e c
x'', p e c
y'')
      | p e c
0 <- p e c
y' = (p e c
x', p e c
y')
      | p e c
1 <- p e c
y' = (p e c
x', p e c
y')
      | Bool
otherwise = (c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
y') p e c
x', c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
y') p e c
y')

-- | Returns an equivalent polynomial from a rational function if possible.
--
-- It returns the numerator if the denominator is a constant.
--
-- >>> toPolynomial $ fromPolynomials (2 * power 1 + 2) (2 :: IndexedPolynomial)
-- Just x + 1
--
-- It returns nothing if the rational function is not equivalent to a polynomial.
--
-- >>> toPolynomial $ fromPolynomials 1 (power 1 :: IndexedPolynomial)
-- Nothing
toPolynomial :: (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) => Function (p e c) -> Maybe (p e c)
toPolynomial :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
Function (p e c) -> Maybe (p e c)
toPolynomial (F p e c
x p e c
1) = p e c -> Maybe (p e c)
forall a. a -> Maybe a
Just p e c
x
toPolynomial (F p e c
x p e c
y)
  | e
0 <- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
y, p e c
y p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0 = p e c -> Maybe (p e c)
forall a. a -> Maybe a
Just (p e c -> Maybe (p e c)) -> p e c -> Maybe (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
y) p e c
x
  | Bool
otherwise = Maybe (p e c)
forall a. Maybe a
Nothing