Skip to content

Commit

Permalink
[Bodigrim#118] Add test to divModE
Browse files Browse the repository at this point in the history
The tests for `divModE` fail, the reason for the function's
incorrect behavior is under investigation.

`divModE'` is a similar division algorith used for
`GaussianIntegers` that works, if no other way to fix the
current process used by `divModE`, `divModE'` will be
used instead.
  • Loading branch information
rockbmb committed Aug 16, 2018
1 parent 3b481b2 commit a4311f3
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 28 deletions.
84 changes: 56 additions & 28 deletions Math/NumberTheory/EisensteinIntegers.hs
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,16 @@ module Math.NumberTheory.EisensteinIntegers
, conjugate
, divE
, divModE
, divModE'
, eisensteinToComplex
, modE
, norm
) where

import qualified Data.Complex as C
import Data.Ratio ((%))
import GHC.Generics (Generic)
import qualified Data.Complex as C
import Data.Ratio ((%))
import Debug.Trace (trace)
import GHC.Generics (Generic)

infix 6 :+

Expand Down Expand Up @@ -79,6 +82,18 @@ instance Num EisensteinInteger where
realEisen :: EisensteinInteger -> Rational
realEisen (a :+ b) = fromInteger a - (b % 2)

divModE'
:: EisensteinInteger
-> EisensteinInteger
-> (EisensteinInteger, EisensteinInteger)
divModE' g h =
let nr :+ ni = g * conjugate h
denom = norm h
q = div nr denom :+ div ni denom
p = h * q
in (q, g - p)


-- | Simultaneous 'div' and 'mod' of Eisenstein integers.
-- The algorithm used here was derived from
-- <http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.892.1219&rep=rep1&type=pdf NTRU over the Eisenstein Integers>
Expand All @@ -87,51 +102,64 @@ divModE
:: EisensteinInteger
-> EisensteinInteger
-> (EisensteinInteger, EisensteinInteger)
divModE alfa beta | norm rho1 < norm rho2 = (gamma1, rho1)
| norm rho1 > norm rho2 = (gamma2, rho2)
divModE alfa beta | norm rho1 < norm rho2 = res1
| norm rho1 > norm rho2 = res2
| norm rho1 == norm rho2 &&
realEisen gamma1 > realEisen gamma2 = (gamma1, rho1)
| otherwise = (gamma2, rho2)
realEisen gamma1 > realEisen gamma2 = res1
| otherwise = res2
where
res1 = trace show' (gamma1, rho1)
res2 = trace show' (gamma2, rho2)
show' = "norm ρ1: " ++ show (norm rho1) ++ "\n" ++
"norm ρ2: " ++ show (norm rho2) ++ "\n" ++
"ρ1': " ++ show rho1' ++ "\n" ++
"ρ2': " ++ show rho2' ++ "\n" ++
"ρ1' * β': " ++ show (beta' * rho1') ++ "\n" ++
"ρ2' * β': " ++ show (beta' * rho2') ++ "\n" ++
"(γ1,ρ1): " ++ show (gamma1, rho1) ++ "\n" ++
"(γ2,ρ2): " ++ show (gamma2, rho2) ++ "\n" ++
"Re(γ1): " ++ show (realEisen gamma1) ++ "\n" ++
"Re(γ2): " ++ show (realEisen gamma2) ++ "\n"

-- @sqrt 3@ is used many times throughout the function, as such it's
-- calculated here once.
sqrt3 :: Double
sqrt3 = sqrt 3

-- First step of assignments performed in the Division Algorithm from
-- Theorem 4.1.1.
alfa' = eisensteinToComplex alfa
beta' = eisensteinToComplex beta
a1 C.:+ b1 = alfa' / beta'
a2 C.:+ b2 = alfa' / beta' - eisensteinToComplex ω
alfa' = eisensteinToComplex alfa
beta' = eisensteinToComplex beta
a1 C.:+ b1 = alfa' / beta'
a2 C.:+ b2 = (alfa' / beta') - eisensteinToComplex ω

-- @round@s a @Double@ and returns that as another @Double@.
dToIToD :: Double -> Double
dToIToD = (fromIntegral :: Integer -> Double) . round
dToIToD :: Double -> Double
dToIToD = (fromIntegral :: Integer -> Double) . round

-- Second step of assignments performed in the Division Algorithm from
-- Theorem 4.1.1.
properFraction' :: Double -> (Integer, Double)
properFraction' = properFraction
rho' :: Double -> Double -> C.Complex Double
rho' a' b' = (snd $ properFraction' a') C.:+ (b' - sqrt3 * dToIToD (b' / sqrt3))
rho1' = rho' a1 b1
rho2' = rho' a2 b2
properFraction' :: Double -> (Integer, Double)
properFraction' = properFraction
rho' :: Double -> Double -> C.Complex Double
rho' a' b' = (snd $ properFraction' a') C.:+ (b' - sqrt3 * dToIToD (b' / sqrt3))
rho1' = rho' a1 b1
rho2' = rho' a2 b2

-- Converts a complex number in the form @a + bι@, where @a, b@ are
-- @Double@s, into an @EisensteinInteger@ in the form @α + βω@, where
-- @α, β@ are @Integer@s.
toEisen :: C.Complex Double -> EisensteinInteger
toEisen (x C.:+ y) = round (x + y / sqrt3) :+ round (y * 2 / sqrt3)
toEisen :: C.Complex Double -> EisensteinInteger
toEisen (x C.:+ y) = round (x + y / sqrt3) :+ round (y * 2 / sqrt3)

-- Third step of assignments performed in the Division Alrgorithm from
-- Third step of assignments performed in the Division Algorithm from
-- Theorem 4.1.1.
rho1 = toEisen $ beta' * rho1'
b1sqrt3' = round $ b1 / sqrt3
gamma1 = ((round a1) + b1sqrt3') :+ (2 * b1sqrt3')
rho2 = toEisen $ beta' * rho2'
b2sqrt3' = round $ b2 / sqrt3
gamma2 = ((round a2) + b2sqrt3') :+ (1 + 2 * b2sqrt3')
rho1 = toEisen $ beta' * rho1'
b1sqrt3' = round $ b1 / sqrt3
gamma1 = ((round a1) + b1sqrt3') :+ (2 * b1sqrt3')
rho2 = toEisen $ beta' * rho2'
b2sqrt3' = round $ b2 / sqrt3
gamma2 = ((round a2) + b2sqrt3') :+ (1 + 2 * b2sqrt3')


-- | Eisenstein integer division, truncating toward negative infinity.
Expand Down
12 changes: 12 additions & 0 deletions test-suite/Math/NumberTheory/EisensteinIntegersTests.hs
Original file line number Diff line number Diff line change
Expand Up @@ -38,9 +38,21 @@ absProperty z = isOrigin || (inFirstSextant && isAssociate)
inFirstSextant = x' > y' && y' >= 0
isAssociate = z' `elem` map (\e -> z * (1 E.:+ 1) ^ e) [0 .. 5]

-- | Verify that `divModE` produces the right quotient and remainder.
divModProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
divModProperty1 x y = (y == 0) || (x `E.divE` y) * y + (x `E.modE` y) == x

-- | Verify that `divModE` produces a remainder smaller than the divisor with
-- regards to the Euclidean domain's function.
modProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
modProperty1 x y = (y == 0) || (E.norm $ x `E.modE` y) < (E.norm y)

testSuite :: TestTree
testSuite = testGroup "EisensteinIntegers" $
[ testSmallAndQuick "forall z . z == signum z * abs z" signumAbsProperty
, testSmallAndQuick "abs z always returns an @EisensteinInteger@ in the\
\ first sextant of the complex plane" absProperty
, testSmallAndQuick "divModE works properly" divModProperty1
, testSmallAndQuick "The remainder's norm is smaller than the divisor's"
modProperty1
]

0 comments on commit a4311f3

Please sign in to comment.