Cbrt.hs

-- |
-- 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 (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
  Zero <> _                    = Zero
  _ <> Zero                    = Zero
  One <> y                     = y
  x <> One                     = x
  Omega <> Omega               = OmegaSquare
  Omega <> OmegaSquare         = One
  OmegaSquare <> Omega         = One
  OmegaSquare <> OmegaSquare   = Omega
  stimes k n = case (k `mod` 3, n) of
    (0, _)           -> One
    (1, symbol)       -> symbol
    (2, Omega)       -> OmegaSquare
    (2, OmegaSquare) -> Omega
    (2, symbol)      -> symbol
    _                -> error "Math.NumberTheory.Moduli.Cbrt: exponentiation undefined."

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

-- | Converts a
-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character cubic symbol>
-- to an Eisenstein Integer.
symbolToNum :: CubicSymbol -> EisensteinInteger
symbolToNum = \case
  Zero        -> 0
  Omega       -> ω
  OmegaSquare -> -1 - ω
  One         -> 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 alpha beta = case beta `A.rem` (1 - ω) of
  -- This checks whether beta is coprime to 3, i.e. divisible by @1 - ω@
  -- In particular, it returns an error if @beta == 0@
  0 -> error "Math.NumberTheory.Moduli.Cbrt: denominator is not coprime to 3."
  _ -> cubicSymbolHelper alpha beta

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

cubicReciprocity :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
-- Note @cubicReciprocity 0 1 = One@. It is better to adopt this convention.
cubicReciprocity _ 1 = One
-- Checks if first argument is zero. Note the second argument is never zero.
cubicReciprocity 0 _ = Zero
-- This checks if the first argument is a unit. Because it's primary,
-- it is enough to pattern match with 1.
cubicReciprocity 1 _ = One
-- Otherwise, cubic reciprocity is called.
cubicReciprocity alpha beta = cubicSymbolHelper beta 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 alpha beta = (gamma, delta, newSymbol)
  where
    newSymbol = stimes (j * m) Omega <> stimes (- m - n) i
    m :+ n = A.quot (delta - 1) 3
    (i, gamma) = getPrimaryDecomposition alphaThreeFree
    (_, delta) = getPrimaryDecomposition beta
    j = wordToInteger jIntWord
    -- This function outputs data such that
    -- @(1 - ω)^jIntWord * alphaThreeFree = alpha@.
    (jIntWord, alphaThreeFree) = splitOff (1 - ω) 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 0 = (Zero, 0)
getPrimaryDecomposition e = case e `A.rem` 3 of
  1            -> (One, e)
  1 :+ 1       -> (OmegaSquare, -ω * e)
  0 :+ 1       -> (Omega, (-1 - ω) * e)
  (-1) :+ 0    -> (One, -e)
  (-1) :+ (-1) -> (OmegaSquare, ω * e)
  0 :+ (-1)    -> (Omega, (1 + ω) * e)
  _            -> error "Math.NumberTheory.Moduli.Cbrt: primary decomposition failed."