fork(5) download
  1. {-# OPTIONS_GHC -O3 -XStrict #-}
  2.  
  3. import Data.Word
  4. import Data.List (sortBy)
  5. import Data.Function (on)
  6.  
  7. nthHam :: Word64 -> (Int, Int, Int)
  8. nthHam n -- n: 1-based 1,2,3...
  9. | n < 2 = case n of
  10. 0 -> error "nthHam: Argument is zero!"
  11. _ -> (0, 0, 0) -- trivial case for 1
  12. | m < 0 = error $ "Not enough triples generated: " ++ show (c,n)
  13. | m >= nb = error $ "Generated band is too narrow: " ++ show (m,nb)
  14. | otherwise = case res of (_, tv) -> tv -- 2^i * 3^j * 5^k
  15. where
  16. lb3 = logBase 2 3; lb5 = logBase 2 5.0
  17. lbrt30 = logBase 2 $ sqrt 30 :: Double -- estimate adjustment as per WP
  18. lg2est = (6 * lb3 * lb5 * fromIntegral n)**(1/3) - lbrt30 -- estimated logval, base 2
  19. (hi,lo) = (lg2est + 1/lg2est, 2 * lg2est - hi) -- hi > log2est > lo
  20. (c, b) = let klmt = floor (hi / lb5)
  21. loopk k ck bndk =
  22. if k > klmt then (ck, bndk) else
  23. let p = hi - fromIntegral k * lb5; jlmt = floor (p / lb3)
  24. loopj j cj bndj =
  25. if j > jlmt then loopk (k + 1) cj bndj else
  26. let q = p - fromIntegral j * lb3
  27. (i, frac) = properFraction q
  28. nj = j + 1; ncj = cj + fromIntegral i + 1
  29. r = hi - frac
  30. nbndj = i `seq` bndj `seq`
  31. if r < lo then bndj
  32. else case (r, (i, j, k)) of
  33. nhd -> nhd `seq` nhd : bndj
  34. in ncj `seq` nbndj `seq` loopj nj ncj nbndj
  35. in loopj 0 ck bndk
  36. in loopk 0 0 []
  37. (m,nb) = ( fromIntegral $ c - n, length b ) -- m 0-based from top, |band|
  38. -- (s,res) = (b, s!!m)
  39. (s,res) = ( sortBy (flip compare `on` fst) b, s!!m ) -- sorted decreasing, result<
  40.  
  41. main = putStrLn $ show $ nthHam 1000000000000
Success #stdin #stdout 0.71s 0KB
stdin
Standard input is empty
stdout
(1126,16930,40)