Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft of CubicSymbol #194

Merged
merged 19 commits into from
May 2, 2020
Merged
Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
27feeae
Draft of CubicSymbol
federico-bongiorno Apr 11, 2020
8c7cb1e
Edited arithmoi.cabal
federico-bongiorno Apr 11, 2020
ad0d206
Extended function to all integers and defined operations for CubicSym…
federico-bongiorno Apr 12, 2020
72ce88f
Addressed comments on pattern matching
federico-bongiorno Apr 12, 2020
d08fc9f
Changed the implementation of the cubicSymbol function and added comm…
federico-bongiorno Apr 13, 2020
d9ea713
Added modular arithmetic, imported functions and error messages.
federico-bongiorno Apr 17, 2020
a05b4b6
Corrected edge cases. Simplified code by adding cubicReciprocity func…
federico-bongiorno Apr 17, 2020
97f4270
Got rid of factoriseBadPrime, changed line spacing and updated comments.
federico-bongiorno Apr 17, 2020
03b7df7
Tidied up code
federico-bongiorno Apr 17, 2020
8ab434b
Improvements to getPrimaryDecomposition function
federico-bongiorno Apr 17, 2020
4730d19
More readable code in getPrimaryDecomposition
federico-bongiorno Apr 17, 2020
eac987b
Added test to check cubicSymbol when the denominator is prime
federico-bongiorno Apr 17, 2020
0533ff6
Minor change to cubicSymbol3 test function
federico-bongiorno Apr 17, 2020
4043c6c
Helper functions return cubic symbols rather than integers
federico-bongiorno Apr 26, 2020
20ab1b8
Added symbolToNum and exponentiation. Changed extractPrimaryContribut…
federico-bongiorno Apr 26, 2020
b00093b
Updating with master
federico-bongiorno Apr 29, 2020
88a218d
Added table lookup for getPrimaryDecomposition
federico-bongiorno Apr 29, 2020
eb9e5ee
Added Haddock comments
federico-bongiorno May 2, 2020
c4ad23b
Changed comments and avoided compiler warning
federico-bongiorno May 2, 2020
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
132 changes: 132 additions & 0 deletions Math/NumberTheory/Moduli/CubicSymbol.hs
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
{-# LANGUAGE LambdaCase #-}

module Math.NumberTheory.Moduli.CubicSymbol
( CubicSymbol(..)
, symbolToNum
, 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.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
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"

symbolToNum :: CubicSymbol -> EisensteinInteger
symbolToNum = \case
Zero -> 0
Omega -> ω
OmegaSquare -> -1 - ω
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 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
federico-bongiorno marked this conversation as resolved.
Show resolved Hide resolved
where
(primaryRemainder, primaryBeta, newSymbol) = extractPrimaryContributions remainder beta
remainder = A.rem alpha beta

-- 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, CubicSymbol)
extractPrimaryContributions alpha beta = (gamma, delta, newSymbol)
where
newSymbol = stimes contribution Omega
contribution = j*m - i*m -i*n
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 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)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you possibly return (CubicSymbol, EisensteinInteger) and make corresponding changes to extractPrimaryContributions and cubicHelper? This will be more natural representation than Integer and Mod 3, I think.

I suggest to pattern match on possible values of remainder and return corresponding cubic symbols directly.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is probably better

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed it but I am disappointed by the result. The code looks messier and it runs more slowly. I thought about it but I could not find any better way of implementing these changes. If you have any suggestions to improve it, let me know.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's surprising, I'll take a look next week. I should revisit M.NT.Q.EisensteinIntegers to check that all operations are as efficient as possible.

In the meantime I have three suggestions:

  • Instead of relying on autoderived stimes (which appears to be awkward) we can actually define our own, supporting zeros and negatives. This basically brings us back to your very first design, sorry for the digression.
instance Semigroup CubicSymbol where
  (<>) = ... 
  stimes k n = case (k `mod` 3, n) of ...
  • Could you please implement symbolToNum :: CubicSymbol -> EisensteinInteger?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did this. In the new commit, the timings went down to the previous level.

-- 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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems we can do this "table lookup" more explicit, pattern-matching by all possible values of remainder. Pattern-matching provides an elegant way to define symb (equal to stimes powerUnit Omega), and then one can choose either e * symbolToNum symb or -e * symbolToNum, depending on which of them is equivalent to 1 modulo 3.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I understand it, this method is not going to work. My understanding is that you want to define factor = e * symbolToNum symb or factor = -e * symbolToNum symb depending on which one is equivalent to 1 mod 3. In general, I don't think that either of them need be equivalent to 1 mod 3. For instance, suppose that I find that (1 + ω) * e A.rem 3 = 1. Then its associated cubic symbol is ω. Then you see that both e * ω and -e * ω cannot have the above property. The issue is that I need to know powerUnit modulo 6 rather that modulo 3. One solution could be to pattern match on remainder and define both factor and symb at the same time. Let me know what you think.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have taken this approach in the latest commit. I would like to make the pattern matching more readable but the program does not compile if I use ω when I list the cases. Also, it gives me a warning for redundant pattern matching which I don't agree with.

inverseRemainder = conjugate remainder
-- Note that this number is the inverse of what is needed.
remainder = e `A.rem` 3
2 changes: 2 additions & 0 deletions arithmoi.cabal
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
75 changes: 75 additions & 0 deletions test-suite/Math/NumberTheory/Moduli/CubicSymbolTests.hs
Original file line number Diff line number Diff line change
@@ -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
]
2 changes: 2 additions & 0 deletions test-suite/Test.hs
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -63,6 +64,7 @@ tests = testGroup "All"
, testGroup "Moduli"
[ ModuliChinese.testSuite
, ModuliClass.testSuite
, ModuliCubic.testSuite
, ModuliDiscreteLogarithm.testSuite
, ModuliEquations.testSuite
, ModuliJacobi.testSuite
Expand Down