{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE LambdaCase #-}
module Math.NumberTheory.Moduli.JacobiSymbol
( JacobiSymbol(..)
, jacobi
, symbolToNum
) where
import Data.Bits
#if __GLASGOW_HASKELL__ < 803
import Data.Semigroup
#endif
import Numeric.Natural
import Math.NumberTheory.Utils
data JacobiSymbol = MinusOne | Zero | One
deriving (JacobiSymbol -> JacobiSymbol -> Bool
(JacobiSymbol -> JacobiSymbol -> Bool)
-> (JacobiSymbol -> JacobiSymbol -> Bool) -> Eq JacobiSymbol
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: JacobiSymbol -> JacobiSymbol -> Bool
$c/= :: JacobiSymbol -> JacobiSymbol -> Bool
== :: JacobiSymbol -> JacobiSymbol -> Bool
$c== :: JacobiSymbol -> JacobiSymbol -> Bool
Eq, Eq JacobiSymbol
Eq JacobiSymbol
-> (JacobiSymbol -> JacobiSymbol -> Ordering)
-> (JacobiSymbol -> JacobiSymbol -> Bool)
-> (JacobiSymbol -> JacobiSymbol -> Bool)
-> (JacobiSymbol -> JacobiSymbol -> Bool)
-> (JacobiSymbol -> JacobiSymbol -> Bool)
-> (JacobiSymbol -> JacobiSymbol -> JacobiSymbol)
-> (JacobiSymbol -> JacobiSymbol -> JacobiSymbol)
-> Ord JacobiSymbol
JacobiSymbol -> JacobiSymbol -> Bool
JacobiSymbol -> JacobiSymbol -> Ordering
JacobiSymbol -> JacobiSymbol -> JacobiSymbol
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: JacobiSymbol -> JacobiSymbol -> JacobiSymbol
$cmin :: JacobiSymbol -> JacobiSymbol -> JacobiSymbol
max :: JacobiSymbol -> JacobiSymbol -> JacobiSymbol
$cmax :: JacobiSymbol -> JacobiSymbol -> JacobiSymbol
>= :: JacobiSymbol -> JacobiSymbol -> Bool
$c>= :: JacobiSymbol -> JacobiSymbol -> Bool
> :: JacobiSymbol -> JacobiSymbol -> Bool
$c> :: JacobiSymbol -> JacobiSymbol -> Bool
<= :: JacobiSymbol -> JacobiSymbol -> Bool
$c<= :: JacobiSymbol -> JacobiSymbol -> Bool
< :: JacobiSymbol -> JacobiSymbol -> Bool
$c< :: JacobiSymbol -> JacobiSymbol -> Bool
compare :: JacobiSymbol -> JacobiSymbol -> Ordering
$ccompare :: JacobiSymbol -> JacobiSymbol -> Ordering
$cp1Ord :: Eq JacobiSymbol
Ord, Int -> JacobiSymbol -> ShowS
[JacobiSymbol] -> ShowS
JacobiSymbol -> String
(Int -> JacobiSymbol -> ShowS)
-> (JacobiSymbol -> String)
-> ([JacobiSymbol] -> ShowS)
-> Show JacobiSymbol
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [JacobiSymbol] -> ShowS
$cshowList :: [JacobiSymbol] -> ShowS
show :: JacobiSymbol -> String
$cshow :: JacobiSymbol -> String
showsPrec :: Int -> JacobiSymbol -> ShowS
$cshowsPrec :: Int -> JacobiSymbol -> ShowS
Show)
instance Semigroup JacobiSymbol where
<> :: JacobiSymbol -> JacobiSymbol -> JacobiSymbol
(<>) = \case
JacobiSymbol
MinusOne -> JacobiSymbol -> JacobiSymbol
negJS
JacobiSymbol
Zero -> JacobiSymbol -> JacobiSymbol -> JacobiSymbol
forall a b. a -> b -> a
const JacobiSymbol
Zero
JacobiSymbol
One -> JacobiSymbol -> JacobiSymbol
forall a. a -> a
id
negJS :: JacobiSymbol -> JacobiSymbol
negJS :: JacobiSymbol -> JacobiSymbol
negJS = \case
JacobiSymbol
MinusOne -> JacobiSymbol
One
JacobiSymbol
Zero -> JacobiSymbol
Zero
JacobiSymbol
One -> JacobiSymbol
MinusOne
{-# SPECIALISE symbolToNum :: JacobiSymbol -> Integer,
JacobiSymbol -> Int,
JacobiSymbol -> Word,
JacobiSymbol -> Natural
#-}
symbolToNum :: Num a => JacobiSymbol -> a
symbolToNum :: JacobiSymbol -> a
symbolToNum = \case
JacobiSymbol
Zero -> a
0
JacobiSymbol
One -> a
1
JacobiSymbol
MinusOne -> -a
1
{-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol,
Natural -> Natural -> JacobiSymbol,
Int -> Int -> JacobiSymbol,
Word -> Word -> JacobiSymbol
#-}
jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi :: a -> a -> JacobiSymbol
jacobi a
_ a
1 = JacobiSymbol
One
jacobi a
a a
b
| a
b a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = String -> JacobiSymbol
forall a. HasCallStack => String -> a
error String
"Math.NumberTheory.Moduli.jacobi: negative denominator"
| a -> Bool
forall a. Bits a => a -> Bool
evenI a
b = String -> JacobiSymbol
forall a. HasCallStack => String -> a
error String
"Math.NumberTheory.Moduli.jacobi: even denominator"
| Bool
otherwise = a -> a -> JacobiSymbol
forall a. (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi' a
a a
b
jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi' :: a -> a -> JacobiSymbol
jacobi' a
0 a
_ = JacobiSymbol
Zero
jacobi' a
1 a
_ = JacobiSymbol
One
jacobi' a
a a
b
| a
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = let n :: JacobiSymbol
n = if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol
MinusOne else JacobiSymbol
One
(Word
z, a
o) = a -> (Word, a)
forall a. Integral a => a -> (Word, a)
shiftToOddCount (a -> a
forall a. Num a => a -> a
negate a
a)
s :: JacobiSymbol
s = if Word -> Bool
forall a. Bits a => a -> Bool
evenI Word
z Bool -> Bool -> Bool
|| a -> Bool
forall a. Bits a => a -> Bool
rem8is1or7 a
b then JacobiSymbol
n else JacobiSymbol -> JacobiSymbol
negJS JacobiSymbol
n
in JacobiSymbol
s JacobiSymbol -> JacobiSymbol -> JacobiSymbol
forall a. Semigroup a => a -> a -> a
<> a -> a -> JacobiSymbol
forall a. (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi' a
o a
b
| a
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
b = case a
a a -> a -> a
forall a. Integral a => a -> a -> a
`rem` a
b of
a
0 -> JacobiSymbol
Zero
a
r -> JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacPS JacobiSymbol
One a
r a
b
| a -> Bool
forall a. Bits a => a -> Bool
evenI a
a = case a -> (Word, a)
forall a. Integral a => a -> (Word, a)
shiftToOddCount a
a of
(Word
z, a
o) -> let r :: JacobiSymbol
r = if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
o Bool -> Bool -> Bool
&& a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol
MinusOne else JacobiSymbol
One
s :: JacobiSymbol
s = if Word -> Bool
forall a. Bits a => a -> Bool
evenI Word
z Bool -> Bool -> Bool
|| a -> Bool
forall a. Bits a => a -> Bool
rem8is1or7 a
b then JacobiSymbol
r else JacobiSymbol -> JacobiSymbol
negJS JacobiSymbol
r
in JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacOL JacobiSymbol
s a
b a
o
| Bool
otherwise = JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacOL (if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
a Bool -> Bool -> Bool
&& a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol
MinusOne else JacobiSymbol
One) a
b a
a
jacPS :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol
jacPS :: JacobiSymbol -> a -> a -> JacobiSymbol
jacPS !JacobiSymbol
acc a
a a
b
| a -> Bool
forall a. Bits a => a -> Bool
evenI a
a = case a -> (Word, a)
forall a. Integral a => a -> (Word, a)
shiftToOddCount a
a of
(Word
z, a
o)
| Word -> Bool
forall a. Bits a => a -> Bool
evenI Word
z Bool -> Bool -> Bool
|| a -> Bool
forall a. Bits a => a -> Bool
rem8is1or7 a
b -> JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacOL (if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
o Bool -> Bool -> Bool
&& a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol -> JacobiSymbol
negJS JacobiSymbol
acc else JacobiSymbol
acc) a
b a
o
| Bool
otherwise -> JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacOL (if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
o Bool -> Bool -> Bool
&& a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol
acc else JacobiSymbol -> JacobiSymbol
negJS JacobiSymbol
acc) a
b a
o
| Bool
otherwise = JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacOL (if a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
a Bool -> Bool -> Bool
&& a -> Bool
forall a. Bits a => a -> Bool
rem4is3 a
b then JacobiSymbol -> JacobiSymbol
negJS JacobiSymbol
acc else JacobiSymbol
acc) a
b a
a
jacOL :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol
jacOL :: JacobiSymbol -> a -> a -> JacobiSymbol
jacOL !JacobiSymbol
acc a
_ a
1 = JacobiSymbol
acc
jacOL !JacobiSymbol
acc a
a a
b = case a
a a -> a -> a
forall a. Integral a => a -> a -> a
`rem` a
b of
a
0 -> JacobiSymbol
Zero
a
r -> JacobiSymbol -> a -> a -> JacobiSymbol
forall a.
(Integral a, Bits a) =>
JacobiSymbol -> a -> a -> JacobiSymbol
jacPS JacobiSymbol
acc a
r a
b
evenI :: Bits a => a -> Bool
evenI :: a -> Bool
evenI a
n = Bool -> Bool
not (a
n a -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
`testBit` Int
0)
rem4is3 :: Bits a => a -> Bool
rem4is3 :: a -> Bool
rem4is3 a
n = a
n a -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
`testBit` Int
1
rem8is1or7 :: Bits a => a -> Bool
rem8is1or7 :: a -> Bool
rem8is1or7 a
n = a
n a -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
`testBit` Int
1 Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== a
n a -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
`testBit` Int
2