Add Artin-Schreier root function.
Add tests.
MCBWM3FBZPTBMBJM4QISAX45C2D75QGHWXPLAYMN3EXJJ2CV4ASAC
L6TQHWB6JUMHD6KZR3RMO23UHQ2E37P2BAFZGDXBL5KWK2WWAIJQC
G25ET6HXYNPM2HJ2L76TGOLVV2LBDA3CIWCAN6SL2WKNTUG56RDQC
TEWWDULPAUYV2VFT3LGUTXGVVDPFFTBRHXQVJAZT63D5Q2NLPS6QC
V7A5CRSQQY3WTZFNH73E6J4YOKE7YJMEW2XW7YJTM37UNTWRMZJQC
5CXVNGYMMIBKPFJFJK7T5HLMUAP2IE7UCFZBFZTY5DBJFRQHPP5QC
VRPOSMITS7VRSIJU6YNEBELCUBTHPMPJD6F6F5PE35R2KECS42KAC
PMF36FUOINXVVNP5XW6ZOFPPWFBRK62O4VCX45NE4Y4GDPL2Z3BQC
6U24OB4HZGUVZZTOEHQKTBOHGV5FKEE4EQDGJF6IP3QQ3HXNNWUQC
P63H7XZA5FGJP5RTASK5NOWQDV4BSKFFCN4AIS7RW44KUK3LK2JQC
KARMFRM7OJF6RDDAVRISRV223RS6AWE6BG33QMFJ4GYLQFZYL56QC
BKOAQCTL55P6HCIDOLA3WKRDJPPZTFDFROKMKRCLXT7ANA3JG4NQC
2Q7SZHYMAFCYLG5MROYQ4BJ7HDIASPLIZJP7SM7BNX4GESUCZXJQC
WFOLCYCZVNRVCEWSDJR36BGSN6IKJ4CMISYCOWBLENMOP2JUTRIQC
O6UAIBEXBNE2XMWQWU6YKGGRLDQFI364JMV6HCYCCI3LBQRD6RPQC
mul (a, b) = mex $ [(nimberProdTable !! a' !! b) `nimberAdd` (nimberProdTable !! a !! b') `nimberAdd` (nimberProdTable !! a' !! b') | a' <- [0 .. a - 1], b' <- [0 .. b - 1]]
mult a b = mex $ [(nimberProdTable !! a' !! b) `nimberAdd` (nimberProdTable !! a !! b') `nimberAdd` (nimberProdTable !! a' !! b') | a' <- [0 .. a - 1], b' <- [0 .. b - 1]]
-- prop_def_add :: Small Int -> Small Int -> Bool
-- prop_def_add Small {getSmall = a} Small {getSmall = b} = abs a `nimberAdd` abs b == fromIntegral (getNimber $ fromIntegral a + fromIntegral b)
prop_def_add :: Bool
prop_def_add = and $ do
i <- [0 .. 15]
j <- [0 .. 15]
pure $ fromIntegral @Int @Nimber (nimberAdd i j) == fromIntegral i + fromIntegral j
-- prop_def_mul :: Small Int -> Small Int -> Bool
-- prop_def_mul Small {getSmall = a} Small {getSmall = b} = abs a `nimberMul` abs b == fromIntegral (getNimber $ fromIntegral a * fromIntegral b)
prop_def_mult :: Bool
prop_def_mult = and $ do
i <- [0 .. 15]
j <- [0 .. 15]
pure $ fromIntegral @Int @Nimber (nimberMult i j) == fromIntegral i * fromIntegral j
semiD = bit (bit m - 1) -- semimultiple of D
a1 = a `shiftR` bit m -- a = a1D+a2
a2 = a .^. (a1 `shiftL` bit m)
b1 = b `shiftR` bit m -- b = b1D+b2
b2 = b .^. (b1 `shiftL` bit m)
c = a2 * b2
in (((a1 + a2) * (b1 + b2) - c) `shiftL` bit m) + a1 * b1 * semiD + c
mult' _ 0 _ = 0
mult' _ _ 0 = 0
mult' _ 1 t = t
mult' _ s 1 = s
mult' k s t =
let semiD = bit (bit k - 1) -- semimultiple of D
s1 = s `shiftR` bit k -- a = a1D+a2
s2 = s .^. (s1 `shiftL` bit k)
t1 = t `shiftR` bit k -- b = b1D+b2
t2 = t .^. (t1 `shiftL` bit k)
c = mult' (k-1) s2 t2
in ((mult' (k-1) (s1 + s2) (t1 + t2) - c) `shiftL` bit k) + mult' (k-1) (mult' (k-1) s1 t1) semiD + c
in mult' m a b
a = n `shiftR` bit m -- n = aD+b
aD = a `shiftL` bit m
b = n .^. aD
semiD = bit (bit m - 1) -- semimultiple of D
sqra = sqr a
in sqra `shiftL` bit m + sqra * semiD + sqr b
sqr' _ 0 = 0
sqr' _ 1 = 1
sqr' k l =
let a = l `shiftR` bit k -- n = aD+b
aD = a `shiftL` bit k
b = l .^. aD
semiD = bit (bit k - 1) -- semimultiple of D
sqra = sqr' (k - 1) a
in sqra `shiftL` bit k + sqra * semiD + sqr' (k - 1) b
in sqr' m n
pow :: (Integral a) => Nimber -> a -> Nimber
_ `pow` 0 = 1
a `pow` n
| n < 0 = recip a `pow` negate n
| otherwise = if even n then sqr a `pow` quot n 2 else powAcc (sqr a) (quot n 2) a
where
powAcc _ 0 x = x
powAcc b m x = if even m then powAcc (sqr b) (quot m 2) x else powAcc (sqr b) (quot m 2) (b * x)
pow :: (Integral a, Bits a) => Nimber -> a -> Nimber
x `pow` n
| n < 0 = recip x `pow` negate n
| otherwise = product . fmap snd . filter (testBit n . fst) . zip [0 ..] . take (1 + floorLog (n + 1)) $ iterate sqr x
a = n `shiftR` bit m -- n = aD+b
aD = a `shiftL` bit m
b = n .^. aD
semiD = bit (bit m - 1) -- semimultiple of D
in (aD + a + b) / (semiD * sqr a + b * (a + b))
recip' _ 1 = 1
recip' k l =
let a = l `shiftR` bit k -- n = aD+b
aD = a `shiftL` bit k
b = l .^. aD
semiD = bit (bit k - 1) -- semimultiple of D
in (aD + a + b) * recip' (k - 1) (semiD * sqr a + b * (a + b))
in recip' m n
a = n `shiftR` bit m -- n = aD+b
aD = a `shiftL` bit m
b = n .^. aD
semiD = bit (bit m - 1) -- semimultiple of D
sqrta = sqrt a
in sqrta `shiftL` bit m + sqrta * sqrt semiD + sqrt b
sqrt' _ 0 = 0
sqrt' _ 1 = 1
sqrt' k l =
let a = l `shiftR` bit k -- n = aD+b
aD = a `shiftL` bit k
b = l .^. aD
semiD = bit (bit k - 1) -- semimultiple of D
sqrta = sqrt' (k - 1) a
in sqrta `shiftL` bit k + sqrta * sqrt semiD + sqrt' (k - 1) b
in sqrt' m n
evens :: [a] -> [a]
evens (x : _ : xs) = x : evens xs
evens xs = xs
-- | @'artinSchreierRoot' n@ is the smallest solution to the equation \(x^2 - x = n\).
-- The algorithm is due to Chin-Long Chen: <https://ieeexplore.ieee.org/document/1056557>.
-- In fields of characteristic 2, the standard quadratic formula does not work, but any quadratic equation can be solved using square roots and Artin-Schreier roots.
--
-- This function is __much__ slower than @'sqrt'@.
artinSchreierRoot :: Nimber -> Nimber
artinSchreierRoot 0 = 0
artinSchreierRoot 1 = 2
artinSchreierRoot 2 = 4
artinSchreierRoot 3 = 6
artinSchreierRoot n =
let m = 1 + floorLog @Int (floorLog n) -- 2^2^m is the order of the smallest field containing n
m' = if n < bit (bit m - 1) then m else m + 1 -- 2^2^m' is the order of the smallest field containing the Artin-Schreier root of n
squares = iterate sqr n
quarts = evens squares
t4k = sum $ take (bit (m' - 1)) quarts -- trace of the Artin-Screier root of n
in if t4k == 1
then
let s = sum $ do
j <- [1 .. bit (m' - 2) - 1]
i <- [j .. bit (m' - 2) - 1]
pure $ (squares !! (shiftL i 1 - 1 .^. bit (m' - 1))) * (squares !! (shiftL j 1 - 2))
in flip clearBit 0 $
s
+ sqr s
+ (squares !! (bit m' - 1))
* (1 + sum (take (bit (m' - 2)) $ drop (bit (m' - 2)) quarts))
else
let y = bit $ bit m' - 1
z = artinSchreierRoot $ sqr y + y + n
in y + z