import Control.Monad
import Data.Bits

import Data.Vector qualified as V
import Data.Vector.Mutable qualified as M

--   Division is implemented using gaussian elimination.

  recip n@Nimber {getNimber = s} =
    let m = 1 + foldl max 0 (S.unions s) -- 2^2^m is the order of the smallest field containing n
        sw i j u v = do
          a <- M.read v i
          b <- M.read v j
          if b `testBit` (2 ^ m - 1 - i)
            then do M.write v i b
                    M.write v j a
                    a' <- M.read u i
                    b' <- M.read u j
                    M.write u i b'
                    M.write u j a'
            else sw i (j + 1) u v
        pivot :: (Bits a, M.PrimMonad m) => M.MVector (M.PrimState m) a -> M.MVector (M.PrimState m) a -> Int -> m ()
        pivot u v i = do sw i i u v
                         p <- M.read u i
                         q <- M.read v i
                         forM_ [i+1..2^m-1] $ \j -> do a <- M.read v j
                                                       when (a `testBit` (2^m-1-i)) $ M.modify v (xor q) j >> M.modify u (xor p) j
        r = V.create $
          do
            tab1 <- M.generate (2 ^ m) (2 ^) -- All two-powers less than 2^2^m
            tab2 <- M.generate (2 ^ m) (nimberToNatural . (n *) . fromIntegral . (2 ^)) -- n times all two-powers less than 2^2^m
            forM_ [0..2^m-1] $ pivot tab1 tab2
            pure tab1
     in fromIntegral $ r V.! (2 ^ m - 1)

  recip 1 = 1
  recip Nimber {getNimber = s} =
    let m     = foldl max 0 $ S.unions s          -- D = 2^2^m is the largest Fermat 2-power less than or equal to n
        aD    = Nimber $ S.filter (S.member m) s  -- n = aD+b
        b     = Nimber $ S.filter (S.notMember m) s
        a     = Nimber $ S.map (S.delete m) $ getNimber aD
        semiD = Nimber . S.singleton $ S.fromList [0..m-1] -- semimultiple of D
     in (aD+a+b) / (semiD * a^2 + a*b + b^2)

                      vector,