    -- |
    -- Module:      Math.NumberTheory.Moduli.Cbrt
    -- Copyright:   (c) 2020 Federico Bongiorno
    -- Licence:     MIT
    -- Maintainer:  Federico Bongiorno <federicobongiorno97@gmail.com>
    --
    -- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character Cubic symbol>
    -- of two Eisenstein Integers.

    {-# LANGUAGE LambdaCase #-}

    module Math.NumberTheory.Moduli.Cbrt
      ( CubicSymbol(..)
      , cubicSymbol
      , symbolToNum
      ) where

    import Math.NumberTheory.Quadratic.EisensteinIntegers
    import Math.NumberTheory.Utils.FromIntegral
    import qualified Data.Euclidean as A
    import Math.NumberTheory.Utils
    import Data.Semigroup

    -- | Represents the
    -- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character cubic residue character>
    -- It is either @0@, @ω@, @ω²@ or @1@.
    data CubicSymbol = Zero | Omega | OmegaSquare | One deriving (CubicSymbol -> CubicSymbol -> Bool
    (CubicSymbol -> CubicSymbol -> Bool)
    -> (CubicSymbol -> CubicSymbol -> Bool) -> Eq CubicSymbol
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: CubicSymbol -> CubicSymbol -> Bool
    == :: CubicSymbol -> CubicSymbol -> Bool
    $c/= :: CubicSymbol -> CubicSymbol -> Bool
    /= :: CubicSymbol -> CubicSymbol -> Bool
    Eq)

    -- | The set of cubic symbols form a semigroup. Note `stimes`
    -- is allowed to take non-positive values. In other words, the set
    -- of non-zero cubic symbols is regarded as a group.
    --
    -- >>> import Data.Semigroup
    -- >>> stimes (-1) Omega
    -- ω²
    -- >>> stimes 0 Zero
    -- 1
    instance Semigroup CubicSymbol where
      CubicSymbol
    Zero <> :: CubicSymbol -> CubicSymbol -> CubicSymbol
    <> CubicSymbol
    _                    = CubicSymbol
    Zero
      CubicSymbol
    _ <> CubicSymbol
    Zero                    = CubicSymbol
    Zero
      CubicSymbol
    One <> CubicSymbol
    y                     = CubicSymbol
    y
      CubicSymbol
    x <> CubicSymbol
    One                     = CubicSymbol
    x
      CubicSymbol
    Omega <> CubicSymbol
    Omega               = CubicSymbol
    OmegaSquare
      CubicSymbol
    Omega <> CubicSymbol
    OmegaSquare         = CubicSymbol
    One
      CubicSymbol
    OmegaSquare <> CubicSymbol
    Omega         = CubicSymbol
    One
      CubicSymbol
    OmegaSquare <> CubicSymbol
    OmegaSquare   = CubicSymbol
    Omega
      stimes :: forall b. Integral b => b -> CubicSymbol -> CubicSymbol
    stimes b
    k CubicSymbol
    n = case (b
    k b -> b -> b
    forall a. Integral a => a -> a -> a
    `mod` b
    3, CubicSymbol
    n) of
        (b
    0, CubicSymbol
    _)           -> CubicSymbol
    One
        (b
    1, CubicSymbol
    symbol)       -> CubicSymbol
    symbol
        (b
    2, CubicSymbol
    Omega)       -> CubicSymbol
    OmegaSquare
        (b
    2, CubicSymbol
    OmegaSquare) -> CubicSymbol
    Omega
        (b
    2, CubicSymbol
    symbol)      -> CubicSymbol
    symbol
        (b, CubicSymbol)
    _                -> [Char] -> CubicSymbol
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "Math.NumberTheory.Moduli.Cbrt: exponentiation undefined."

    instance Show CubicSymbol where
      show :: CubicSymbol -> [Char]
    show = \case
        CubicSymbol
    Zero         -> [Char]
    "0"
        CubicSymbol
    Omega        -> [Char]
    "ω"
        CubicSymbol
    OmegaSquare  -> [Char]
    "ω²"
        CubicSymbol
    One          -> [Char]
    "1"

    -- | Converts a
    -- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character cubic symbol>
    -- to an Eisenstein Integer.
    symbolToNum :: CubicSymbol -> EisensteinInteger
    symbolToNum :: CubicSymbol -> EisensteinInteger
    symbolToNum = \case
      CubicSymbol
    Zero        -> EisensteinInteger
    0
      CubicSymbol
    Omega       -> EisensteinInteger
    ω
      CubicSymbol
    OmegaSquare -> -EisensteinInteger
    1 EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    - EisensteinInteger
    ω
      CubicSymbol
    One         -> EisensteinInteger
    1

    -- The algorithm `cubicSymbol` is adapted from
    -- <https://cs.au.dk/~gudmund/Documents/cubicres.pdf here>.
    -- It is divided in the following steps.
    --
    -- (1) Check whether @beta@ is coprime to 3.
    -- (2) Replace @alpha@ by the remainder of @alpha@ mod @beta@
    --     This does not affect the cubic symbol.
    -- (3) Replace @alpha@ and @beta@ by their associated primary
    --     divisors and keep track of how their cubic residue changes.
    -- (4) Check if any of the two numbers is a zero or a unit. In this
    --     case, return their cubic residue.
    -- (5) Otherwise, invoke cubic reciprocity by swapping @alpha@ and
    --     @beta@. Note both numbers have to be primary.
    --     Return to Step 2.

    -- | <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character Cubic symbol>
    -- of two Eisenstein Integers.
    -- The first argument is the numerator and the second argument
    -- is the denominator. The latter must be coprime to @3@.
    -- This condition is checked.
    --
    -- If the arguments have a common factor, the result
    -- is 'Zero', otherwise it is either 'Omega', 'OmegaSquare' or 'One'.
    --
    -- >>> cubicSymbol (45 + 23*ω) (11 - 30*ω)
    -- 0
    -- >>> cubicSymbol (31 - ω) (1 +10*ω)
    -- ω
    cubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbol EisensteinInteger
    alpha EisensteinInteger
    beta = case EisensteinInteger
    beta EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Euclidean a => a -> a -> a
    `A.rem` (EisensteinInteger
    1 EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    - EisensteinInteger
    ω) of
      -- This checks whether beta is coprime to 3, i.e. divisible by @1 - ω@
      -- In particular, it returns an error if @beta == 0@
      EisensteinInteger
    0 -> [Char] -> CubicSymbol
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "Math.NumberTheory.Moduli.Cbrt: denominator is not coprime to 3."
      EisensteinInteger
    _ -> EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbolHelper EisensteinInteger
    alpha EisensteinInteger
    beta

    cubicSymbolHelper :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbolHelper :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbolHelper EisensteinInteger
    alpha EisensteinInteger
    beta = EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicReciprocity EisensteinInteger
    primaryRemainder EisensteinInteger
    primaryBeta CubicSymbol -> CubicSymbol -> CubicSymbol
    forall a. Semigroup a => a -> a -> a
    <> CubicSymbol
    newSymbol
      where
        (EisensteinInteger
    primaryRemainder, EisensteinInteger
    primaryBeta, CubicSymbol
    newSymbol) = EisensteinInteger
    -> EisensteinInteger
    -> (EisensteinInteger, EisensteinInteger, CubicSymbol)
    extractPrimaryContributions EisensteinInteger
    remainder EisensteinInteger
    beta
        remainder :: EisensteinInteger
    remainder = EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Euclidean a => a -> a -> a
    A.rem EisensteinInteger
    alpha EisensteinInteger
    beta

    cubicReciprocity :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    -- Note @cubicReciprocity 0 1 = One@. It is better to adopt this convention.
    cubicReciprocity :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicReciprocity EisensteinInteger
    _ EisensteinInteger
    1 = CubicSymbol
    One
    -- Checks if first argument is zero. Note the second argument is never zero.
    cubicReciprocity EisensteinInteger
    0 EisensteinInteger
    _ = CubicSymbol
    Zero
    -- This checks if the first argument is a unit. Because it's primary,
    -- it is enough to pattern match with 1.
    cubicReciprocity EisensteinInteger
    1 EisensteinInteger
    _ = CubicSymbol
    One
    -- Otherwise, cubic reciprocity is called.
    cubicReciprocity EisensteinInteger
    alpha EisensteinInteger
    beta = EisensteinInteger -> EisensteinInteger -> CubicSymbol
    cubicSymbolHelper EisensteinInteger
    beta EisensteinInteger
    alpha

    -- | This function takes two Eisenstein intgers @alpha@ and @beta@ and returns
    -- three arguments @(gamma, delta, newSymbol)@. @gamma@ and @delta@ are the
    -- associated primary numbers of alpha and beta respectively. @newSymbol@
    -- is the cubic symbol measuring the discrepancy between the cubic residue
    -- of @alpha@ and @beta@, and the cubic residue of @gamma@ and @delta@.
    extractPrimaryContributions :: EisensteinInteger -> EisensteinInteger -> (EisensteinInteger, EisensteinInteger, CubicSymbol)
    extractPrimaryContributions :: EisensteinInteger
    -> EisensteinInteger
    -> (EisensteinInteger, EisensteinInteger, CubicSymbol)
    extractPrimaryContributions EisensteinInteger
    alpha EisensteinInteger
    beta = (EisensteinInteger
    gamma, EisensteinInteger
    delta, CubicSymbol
    newSymbol)
      where
        newSymbol :: CubicSymbol
    newSymbol = Integer -> CubicSymbol -> CubicSymbol
    forall b. Integral b => b -> CubicSymbol -> CubicSymbol
    forall a b. (Semigroup a, Integral b) => b -> a -> a
    stimes (Integer
    j Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer
    m) CubicSymbol
    Omega CubicSymbol -> CubicSymbol -> CubicSymbol
    forall a. Semigroup a => a -> a -> a
    <> Integer -> CubicSymbol -> CubicSymbol
    forall b. Integral b => b -> CubicSymbol -> CubicSymbol
    forall a b. (Semigroup a, Integral b) => b -> a -> a
    stimes (- Integer
    m Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    - Integer
    n) CubicSymbol
    i
        Integer
    m :+ Integer
    n = EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Euclidean a => a -> a -> a
    A.quot (EisensteinInteger
    delta EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    - EisensteinInteger
    1) EisensteinInteger
    3
        (CubicSymbol
    i, EisensteinInteger
    gamma) = EisensteinInteger -> (CubicSymbol, EisensteinInteger)
    getPrimaryDecomposition EisensteinInteger
    alphaThreeFree
        (CubicSymbol
    _, EisensteinInteger
    delta) = EisensteinInteger -> (CubicSymbol, EisensteinInteger)
    getPrimaryDecomposition EisensteinInteger
    beta
        j :: Integer
    j = Word -> Integer
    wordToInteger Word
    jIntWord
        -- This function outputs data such that
        -- @(1 - ω)^jIntWord * alphaThreeFree = alpha@.
        (Word
    jIntWord, EisensteinInteger
    alphaThreeFree) = EisensteinInteger -> EisensteinInteger -> (Word, EisensteinInteger)
    forall a. (Eq a, GcdDomain a) => a -> a -> (Word, a)
    splitOff (EisensteinInteger
    1 EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    - EisensteinInteger
    ω) EisensteinInteger
    alpha

    -- | This function takes an Eisenstein number @e@ and returns @(symbol, delta)@
    -- where @delta@ is its associated primary integer and @symbol@ is the
    -- cubic symbol discrepancy between @e@ and @delta@. @delta@ is defined to be
    -- the unique associated Eisenstein Integer to @e@ such that
    -- \( \textrm{delta} \equiv 1 (\textrm{mod} 3) \).
    -- Note that @delta@ exists if and only if @e@ is coprime to 3. In this
    -- case, an error message is displayed.
    getPrimaryDecomposition :: EisensteinInteger -> (CubicSymbol, EisensteinInteger)
    -- This is the case where a common factor between @alpha@ and @beta@ is detected.
    -- In this instance @cubicReciprocity@ will return `Zero`.
    -- Strictly speaking, this is not a primary decomposition.
    getPrimaryDecomposition :: EisensteinInteger -> (CubicSymbol, EisensteinInteger)
    getPrimaryDecomposition EisensteinInteger
    0 = (CubicSymbol
    Zero, EisensteinInteger
    0)
    getPrimaryDecomposition EisensteinInteger
    e = case EisensteinInteger
    e EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Euclidean a => a -> a -> a
    `A.rem` EisensteinInteger
    3 of
      EisensteinInteger
    1            -> (CubicSymbol
    One, EisensteinInteger
    e)
      Integer
    1 :+ Integer
    1       -> (CubicSymbol
    OmegaSquare, -EisensteinInteger
    ω EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    * EisensteinInteger
    e)
      Integer
    0 :+ Integer
    1       -> (CubicSymbol
    Omega, (-EisensteinInteger
    1 EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    - EisensteinInteger
    ω) EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    * EisensteinInteger
    e)
      (-1) :+ Integer
    0    -> (CubicSymbol
    One, -EisensteinInteger
    e)
      (-1) :+ (-1) -> (CubicSymbol
    OmegaSquare, EisensteinInteger
    ω EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    * EisensteinInteger
    e)
      Integer
    0 :+ (-1)    -> (CubicSymbol
    Omega, (EisensteinInteger
    1 EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    + EisensteinInteger
    ω) EisensteinInteger -> EisensteinInteger -> EisensteinInteger
    forall a. Num a => a -> a -> a
    * EisensteinInteger
    e)
      EisensteinInteger
    _            -> [Char] -> (CubicSymbol, EisensteinInteger)
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "Math.NumberTheory.Moduli.Cbrt: primary decomposition failed."