Draft of CubicSymbol

Math/NumberTheory/Moduli/CubicSymbol.hs

@@ -0,0 +1,125 @@
+{-# LANGUAGE LambdaCase #-}
+{-# LANGUAGE DataKinds #-}
+
+module Math.NumberTheory.Moduli.CubicSymbol
+  ( CubicSymbol(..)
+  , cubicSymbol
+  ) where
+
+import Math.NumberTheory.Quadratic.EisensteinIntegers
+import Math.NumberTheory.Utils.FromIntegral
+import qualified Data.Euclidean as A
+import Math.NumberTheory.Utils
+import Data.Semigroup
+import Data.Mod.Word
+import Data.Maybe
+import Data.List
+
+data CubicSymbol = Zero | Omega | OmegaSquare | One deriving (Eq)
+
+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
+
+instance Show CubicSymbol where
+  show = \case
+    Zero         -> "0"
+    Omega        -> "ω"
+    OmegaSquare  -> "ω²"
+    One          -> "1"
+
+-- The algorithm cubicSymbol takes two Eisentein numbers @alpha@ and @beta@ and returns
+-- their cubic residue. 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. If it
+--    is, return their cubic residue.
+-- 5) If not, invoke cubic reciprocity by swapping @alpha@ and
+--    @beta@. Note both numbers have to be primary.
+--    Return to Step 2.
+
+-- This function takes two Eisenstein integers and returns their cubic residue character.
+-- Note that the second argument must be coprime to 3 else the algorithm returns an error.
+cubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
+cubicSymbol alpha beta = case betaNorm `mod` 3 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
+  where
+    betaNorm = norm beta
+
+cubicSymbolHelper :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
+cubicSymbolHelper alpha beta = cubicReciprocity primaryRemainder primaryBeta <> newSymbol
+  where
+    (primaryRemainder, primaryBeta, symbolExponent) = extractPrimaryContributions remainder beta
+    remainder = A.rem alpha beta
+    newSymbol = exponentiation unmodularExponent Omega
+    unmodularExponent = wordToInt (unMod symbolExponent)
+    exponentiation k x = if k == 0 then One else stimes k x
+
+-- This function first checks if its arguments are zeros or units. If they are not,
+-- it invokes cubic reciprocity by calling cubicSymbolHelper with swapped arguments.
+cubicReciprocity :: EisensteinInteger -> EisensteinInteger -> CubicSymbol
+-- Note @cubicReciprocity 0 1 = One@. It turns out 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, contribution)@. @gamma@ and @delta@ are the associated
+-- primary numbers to alpha and beta respectively. @contribution@ is a an integer
+-- defined mod 3 which measures the difference between the cubic residue of @alpha@
+-- and @beta@ with respect to the cubic residue of @gamma@ and @delta@.
+extractPrimaryContributions :: EisensteinInteger -> EisensteinInteger -> (EisensteinInteger, EisensteinInteger, Mod 3)
+extractPrimaryContributions alpha beta = (gamma, delta, contribution)
+  where
+    contribution = j*m - i*m -i*n
+    [i, j, m, n] = map fromIntegral [iInt, jInt, mInt, nInt]
+    mInt :+ nInt = A.quot (delta - 1) 3
+    (iInt, gamma) = getPrimaryDecomposition alphaThreeFree
+    (_, delta) = getPrimaryDecomposition beta
+    jInt = wordToInteger jIntWord
+    -- This function outputs data such that
+    -- @(1 - ω)^jIntWord * alphaThreeFree = alpha@.
+    (jIntWord, alphaThreeFree) = splitOff (1 - ω) alpha
+
+-- This function takes an Eisenstein number and returns its primary decomposition @(powerUnit, factor)@
+-- That is, given @e@ coprime with 3, it returns a unique integer (mod 6) @powerUnit@ and a unique
+-- Eisenstein number @factor@ such that @(1 + ω)^powerUnit * e = 1 + 3*factor@.
+-- Note that L.findIndex cannot return Nothing. This happens only if @e@ is not
+-- coprime with 3. This cannot happen since @U.splitOff@ is called just before.
+getPrimaryDecomposition :: EisensteinInteger -> (Integer, 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 = (0, 0)
+getPrimaryDecomposition e = (toInteger powerUnit, factor)
+  where
+    factor = unit * e
+    unit = (1 + ω)^powerUnit
+    powerUnit = fromMaybe
+      (error "Math.NumberTheory.Moduli.CubicSymbol: primary decomposition failed.")
+      findPowerUnit
+    -- Note that the units in @ids@ are ordered in the following way:
+    -- The i^th element of @ids@ is @(1 + ω)^i@ starting from i = 0@
+    -- That is the i^th unit counting anticlockwise starting with 1.
+    findPowerUnit = elemIndex inverseRemainder ids
+    inverseRemainder = conjugate remainder
+    -- Note that this number is the inverse of what is needed.
+    remainder = e `A.rem` 3

arithmoi.cabal

@@ -59,6 +59,7 @@ library
     Math.NumberTheory.Moduli
     Math.NumberTheory.Moduli.Chinese
     Math.NumberTheory.Moduli.Class
+    Math.NumberTheory.Moduli.CubicSymbol
     Math.NumberTheory.Moduli.DiscreteLogarithm
     Math.NumberTheory.Moduli.Equations
     Math.NumberTheory.Moduli.Jacobi
@@ -147,6 +148,7 @@ test-suite spec
     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,75 @@
+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 the numerators
+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 the denominators
+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 the cubic symbol is correct when the denominator is primebeta
+-- as explanined in § 3.3.2 in https://en.wikipedia.org/wiki/Cubic_reciprocity
+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