Compute BCH-bounds for all cyclic codes of length 31 over F2

orbits.hs

import Data.List (sort, nub)

type Power = Int
type IndexCode = Int
type Flag = Int
type Orbit = [Power]
type Roots = [Power]

type P = Int
type N = Int

-- Next item in a cyclotomic orbit of given element.
next :: N -> P -> Power -> Power
next n p pow = p*pow `mod` n

orbit :: N -> P -> Power -> Orbit
orbit n p start = orbit' n p start (next n p start) []
  where
    orbit' n p start current list = if start == current
      then current:list
      else orbit' n p start (next n p current) (current:list)

initialOrbits :: [Orbit]
initialOrbits = [0] : (nub $ map (sort . orbit 31 2) [1..30])
-- initialOrbits = [ [0]
--                 , [1,2,4,8,16]
--                 , [3,6,12,17,24]
--                 , [5,9,10,18,20]
--                 , [7,14,19,25,28]
--                 , [11,13,21,22,26]
--                 , [15,23,27,29,30]        
--                 ]

-- Two functions to decode binary characteristic set of subset of indices from numeric representation.
gorner :: Power -> [Power] -> [Flag] -> [Power]
gorner _ p [] = p
gorner current accum (f:fs) = if f == 1
  then gorner (current + 1) (current:accum) fs
  else gorner (current + 1) accum fs

extract :: IndexCode -> [Power]
extract n = gorner 0 [] $ extract' n []
  where
    extract' n accum = if n == 0
      then accum
      else extract' (n `div` 2) (accum ++ [n `mod` 2])

-- Glue together set of roots of a divisor
buildRootsSet :: [Orbit] -> IndexCode -> Roots
buildRootsSet selection idx = sort $ concatMap ((!!) selection) (extract idx)

-- Compute BCH-bound
-- Case of the code with no nonzero words is not considered as special here.
-- In a result there will be single code with bound 32. It is it.
-- It is evident since it corresponds to x^{31}-1 which generates zero ideal in a ring we work in -- F[x]/(x^{31}-1). 
countEstimate :: Roots -> Int
countEstimate rs = countEstimate' 0 0 (-1) rs + 1
  where
    countEstimate' m accum last [] = max m accum
    countEstimate' m accum last (r:rs) = if r - last == 1
      then countEstimate' (max m accum) (accum + 1) r rs
      else countEstimate' (max m accum) 1 r rs

main :: IO ()
main = print results
  where estimates = map (countEstimate . buildRootsSet initialOrbits) [0..127]
        results = map (\v -> (v, length (filter (==v) estimates))) (nub estimates)