Judy

judy.hs

import           Control.Monad            (when, forM_)
import qualified Data.Judy                as J
import           Data.Maybe               (fromJust)
import qualified Data.Vector.Unboxed      as V
import           Control.Applicative

p10 n = do
    let r  = floor (sqrt (fromIntegral n))
    let n' = n `div` r - 1
    let v1 = V.generate r  $ \i -> n `div` (i + 1)  -- [n `div` i | i <- [1..r]]
    let v2 = V.generate n' $ \i -> n' - i           -- [n', n' - 1 .. 1]
    let vs = v1 V.++ v2
    s <- J.new :: IO (J.JudyL Int)
    forM_ (V.toList vs) $ \v ->
        J.insert (fromIntegral v) (v * (v + 1) `div` 2 - 1) s
    let loop p
          | p >= r    = pure ()
          | otherwise = do
              Just cur <- J.lookup (fromIntegral p) s
              Just sp  <- J.lookup (fromIntegral p - 1) s
              when (cur > sp) $ do
                  let p2 = p * p
                  V.forM_ (V.takeWhile (>= p2) vs) $ \v -> do
                     let vp = v `div` p
                     Just svp <- J.lookup (fromIntegral vp) s
                     J.adjust (\v -> v - p * (svp - sp)) (fromIntegral v) s
              loop (p + 1)
    loop 2
    fromJust <$> J.lookup (fromIntegral n) s

main = do 
    n <- p10 (2*10^9 :: Int)
    print n