File: Cbrt.hs

package info (click to toggle)
haskell-arithmoi 0.13.2.0-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 964 kB
  • sloc: haskell: 10,379; makefile: 3
file content (160 lines) | stat: -rw-r--r-- 6,649 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
-- |
-- 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."