import Control.Monad
import Control.Monad.ST
import Data.Array.ST
import Data.Maybe

data Primality
	= Prime
	| OneFactor Integer
	| TwoFactor Integer Integer
	| ManyFactor
	deriving (Eq, Ord, Read, Show)

addFactor :: Integer -> Primality -> Primality
addFactor f Prime = OneFactor f
addFactor f (OneFactor f') = TwoFactor f f'
addFactor _ _ = ManyFactor

isSemiprime :: Integer -> Primality -> Bool
isSemiprime n (OneFactor f   ) = f * f  == n
isSemiprime n (TwoFactor f f') = f * f' == n
isSemiprime _ _ = False

nextSemiprime :: Integer -> Integer
nextSemiprime n = runST $ do
	arr <- newSTArray (2,2*n) Prime
	forM_ [2..n] $ \p -> do
		primality <- readArray arr p
		when (primality == Prime) $
			forM_ [2*p,3*p..2*n] $ \i ->
				modifyArray arr i (addFactor p)
	-- going up to 2*n is big enough by Bertrand's postulate
	fromJust <$> findM (\i -> isSemiprime i <$> readArray arr i) [n+1..2*n]

findM :: Monad m => (a -> m Bool) -> [a] -> m (Maybe a)
findM f [] = return Nothing
findM f (x:xs) = do
	done <- f x
	if done then return (Just x) else findM f xs

modifyArray :: (MArray a e m, Ix i) => a i e -> i -> (e -> e) -> m ()
modifyArray arr ix f = readArray arr ix >>= writeArray arr ix . f

-- just to fix the type
newSTArray :: Ix i => (i,i) -> e -> ST s (STArray s i e)
newSTArray = newArray

main = print . take 100 . iterate nextSemiprime $ 4