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