Merge pull request #194 from folidota/folidota/cubic_symbols

CubicSymbol

Math/NumberTheory/Moduli/CubicSymbol.hs

@@ -0,0 +1,159 @@
+-- |
+-- Module:      Math.NumberTheory.Moduli.CubicSymbol
+-- Copyright:   (c) 2020 Federico Bongiorno
+-- Licence:     MIT
+-- Maintainer:  Federico Bongiorno <[email protected]>
+--
+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character Cubic symbol>
+-- of two Eisenstein Integers.
+
+{-# LANGUAGE LambdaCase #-}
+
+module Math.NumberTheory.Moduli.CubicSymbol
+  ( 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.
+--
+-- >>> stimes -1 ω
+-- ω²
+-- >>> stimes 0 0
+-- 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.CubicSymbol: 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.CubicSymbol: 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.CubicSymbol: primary decomposition failed."

arithmoi.cabal

@@ -56,6 +56,7 @@ library
     Math.NumberTheory.Moduli
     Math.NumberTheory.Moduli.Chinese
     Math.NumberTheory.Moduli.Class
+    Math.NumberTheory.Moduli.CubicSymbol
     Math.NumberTheory.Moduli.Equations
     Math.NumberTheory.Moduli.Multiplicative
     Math.NumberTheory.Moduli.Singleton
@@ -135,6 +136,7 @@ test-suite arithmoi-tests
     Math.NumberTheory.Moduli.ChineseTests
     Math.NumberTheory.Moduli.DiscreteLogarithmTests
     Math.NumberTheory.Moduli.ClassTests
+    Math.NumberTheory.Moduli.CubicSymbolTests
     Math.NumberTheory.Moduli.EquationsTests
     Math.NumberTheory.Moduli.JacobiTests
     Math.NumberTheory.Moduli.PrimitiveRootTests

test-suite/Math/NumberTheory/Moduli/CubicSymbolTests.hs

@@ -0,0 +1,89 @@
+-- |
+-- Module:      Math.NumberTheory.Moduli.CubicSymbol
+-- Copyright:   (c) 2020 Federico Bongiorno
+-- Licence:     MIT
+-- Maintainer:  Federico Bongiorno <[email protected]>
+--
+-- Test for Math.NumberTheory.Moduli.CubicSymbol
+--
+
+module Math.NumberTheory.Moduli.CubicSymbolTests
+  ( testSuite
+  ) where
+
+import Math.NumberTheory.Moduli.CubicSymbol
+import Math.NumberTheory.Quadratic.EisensteinIntegers
+import Math.NumberTheory.Primes
+import qualified Data.Euclidean as A
+import Data.List
+import Test.Tasty
+import Math.NumberTheory.TestUtils
+
+-- Checks multiplicative property of numerators. In details,
+-- @cubicSymbol1 alpha1 alpha2 beta@ checks that
+-- @(cubicSymbol alpha1 beta) <> (cubicSymbol alpha2 beta) == (cubicSymbol alpha1*alpha2 beta)@
+cubicSymbol1 :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool
+cubicSymbol1 alpha1 alpha2 beta = isBadDenominator beta || cubicSymbolNumerator alpha1 alpha2 beta
+
+cubicSymbolNumerator :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool
+cubicSymbolNumerator alpha1 alpha2 beta = (symbol1 <> symbol2) == symbolProduct
+  where
+    symbol1 = cubicSymbol alpha1 beta
+    symbol2 = cubicSymbol alpha2 beta
+    symbolProduct = cubicSymbol alphaProduct beta
+    alphaProduct = alpha1 * alpha2
+
+-- Checks multiplicative property of denominators. In details,
+-- @cubicSymbol2 alpha beta1 beta2@ checks that
+-- @(cubicSymbol alpha beta1) <> (cubicSymbol alpha beta2) == (cubicSymbol alpha beta1*beta2)@
+cubicSymbol2 :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool
+cubicSymbol2 alpha beta1 beta2 = isBadDenominator beta1 || isBadDenominator beta2 || cubicSymbolDenominator alpha beta1 beta2
+
+cubicSymbolDenominator :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool
+cubicSymbolDenominator alpha beta1 beta2 = (symbol1 <> symbol2) == symbolProduct
+  where
+    symbol1 = cubicSymbol alpha beta1
+    symbol2 = cubicSymbol alpha beta2
+    symbolProduct = cubicSymbol alpha betaProduct
+    betaProduct = beta1 * beta2
+
+-- Checks that `cubicSymbol` agrees with the computational definition
+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Definition here>
+-- when the denominator is prime.
+cubicSymbol3 :: EisensteinInteger -> Prime EisensteinInteger -> Bool
+cubicSymbol3 alpha prime = isBadDenominator beta || cubicSymbol alpha beta == cubicSymbolPrime alpha beta
+    where beta = unPrime prime
+
+cubicSymbolPrime :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
+cubicSymbolPrime alpha beta = findCubicSymbol residue beta
+  where
+    residue = foldr f 1 listOfAlphas
+    f x y = (x * y) `A.rem` beta
+    listOfAlphas = genericReplicate alphaExponent alpha
+    -- Exponent is defined to be 1/3*(@betaNorm@ - 1).
+    alphaExponent = betaNorm `div` 3
+    betaNorm = norm beta
+
+isBadDenominator :: EisensteinInteger -> Bool
+isBadDenominator x = modularNorm == 0
+  where
+    modularNorm = norm x `mod` 3
+
+-- This complication is necessary because it may happen that the residue field
+-- of @beta@ has characteristic two. In this case 1=-1 and the Euclidean algorithm
+-- can return both. Therefore it is not enough to pattern match for the values
+-- which give a well defined @cubicSymbol@.
+findCubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
+findCubicSymbol residue beta
+  | residue `A.rem` beta == 0             = Zero
+  | (residue - ω) `A.rem` beta == 0       = Omega
+  | (residue + 1 + ω) `A.rem` beta == 0   = OmegaSquare
+  | (residue - 1) `A.rem` beta == 0       = One
+  | otherwise                             = error "Math.NumberTheory.Moduli.CubicSymbol: invalid EisensteinInteger."
+
+testSuite :: TestTree
+testSuite = testGroup "CubicSymbol"
+  [ testSmallAndQuick "multiplicative property of numerators" cubicSymbol1
+  , testSmallAndQuick "multiplicative property of denominators" cubicSymbol2
+  , testSmallAndQuick "cubic residue with prime denominator" cubicSymbol3
+  ]

test-suite/Test.hs

@@ -9,6 +9,7 @@ import qualified Math.NumberTheory.Recurrences.LinearTests as RecurrencesLinear
 
 import qualified Math.NumberTheory.Moduli.ChineseTests as ModuliChinese
 import qualified Math.NumberTheory.Moduli.ClassTests as ModuliClass
+import qualified Math.NumberTheory.Moduli.CubicSymbolTests as ModuliCubic
 import qualified Math.NumberTheory.Moduli.DiscreteLogarithmTests as ModuliDiscreteLogarithm
 import qualified Math.NumberTheory.Moduli.EquationsTests as ModuliEquations
 import qualified Math.NumberTheory.Moduli.JacobiTests as ModuliJacobi
@@ -63,6 +64,7 @@ tests = testGroup "All"
   , testGroup "Moduli"
     [ ModuliChinese.testSuite
     , ModuliClass.testSuite
+    , ModuliCubic.testSuite
     , ModuliDiscreteLogarithm.testSuite
     , ModuliEquations.testSuite
     , ModuliJacobi.testSuite