fork download
  1. {-# OPTIONS -O2 #-}
  2.  
  3. -- directly find n-th Hamming number, in ~ O(n^{2/3}) time
  4. -- based on "top band" idea by Louis Klauder, from DDJ discussion
  5. -- by Will Ness, original post: drdobbs.com/blogs/architecture-and-design/228700538
  6.  
  7. import Data.List -- cf. ideone.com/q3fma
  8. import Data.Function
  9.  
  10. main = let (r,t) = nthHam 1000000000000 in print t -- >> print (trival t)
  11.  
  12. lb3 = logBase 2 3; lb5 = logBase 2 5; lb30_2 = logBase 2 30 / 2
  13. trival (i,j,k) = 2^i * 3^j * 5^k
  14. estval n
  15. | n > 500000 = (v - lb30_2 + (3/v), 6/v) -- the space tweak! (thx, GBG!)
  16. | n > 500000 = (v - 2.4496 , 0.0076 ) -- empirical estimation
  17. | n > 50000 = (v - 2.4424 , 0.0146 ) -- correction, base 2
  18. | n > 500 = (v - 2.3948 , 0.0723 ) -- (dist,width)
  19. | n > 1 = (v - 2.2506 , 0.2887 ) -- around (log $ sqrt 30),
  20. | otherwise = (v - 2.2506 , 0.5771 ) -- says WP
  21. where v = (6*lb3*lb5* fromIntegral n)**(1/3) -- estimated logval, base 2
  22.  
  23. nthHam :: Int -> (Double, (Int, Int, Int)) -- ( 64bit: use Int!!! NB! )
  24. nthHam n -- n: 1-based: 1,2,3...
  25. | n <= 0 = error $ "n is 1--based: must be n > 0: " ++ show n
  26. | w >= 1 = error $ "Breach of contract: (w < 1): " ++ show w
  27. | m < 0 = error $ "Not enough triples generated: " ++ show (c,n)
  28. | m >= nb = error $ "Generated band is too narrow: " ++ show (m,nb)
  29. | otherwise = sortBy (flip compare `on` fst) b !! m -- m-th from top in sorted band
  30. where
  31. (hi,w) = estval n -- hi > logval > hi-w
  32. m = fromIntegral (c - n) -- target index, from top
  33. nb = length b -- length of the band
  34. {- (Sum c,b) = fold [ (Sum (fromIntegral i+1), ... ] -}
  35. (c,b) = foldl_ (\(c,b) (i,t)-> let c2=c+i in c2`seq` -- ( total count, the band )
  36. case t of []-> (c2,b);[v]->(c2,v:b) ) (0,[])
  37. -- prod (sum,concat) . unzip $
  38. [ ( fromIntegral i+1, -- total triples w/ this (j,k)
  39. [ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band
  40. | k <- [ 0 .. floor ( hi /lb5) ], let p = fromIntegral k*lb5,
  41. j <- [ 0 .. floor ((hi-p)/lb3) ], let q = fromIntegral j*lb3 + p,
  42. let (i,frac) = pr (hi-q) ; r = hi - frac -- r = i + q
  43. ] where pr = properFraction -- pr 1.24 => (1,0.24)
  44. prod (f,g) (x,y) = (f x,g y) -- (f *** g)
  45. foldl_ = foldl'
Success #stdin #stdout 0.7s 11820KB
stdin
 
2020-03-23: 
now using Int in nthHam :: Int -> (Double, (Int, Int, Int)) 
instead of the previous Integer:
    nthHam :: Integer -> (Double, (Int, Int, Int)) 

                                      Oct-2017: 
1T:              (1126,16930,40)     1.79s    n^0.67
100B:            (10178,1384,279)    0.38s
20B: 0.80s-5.8MB (4077,927,890)      0.14s
10B: 0.45s-5.8MB (2177,8,1659)       0.08s
1B:  0.05s-4.8MB (1334,335,404)
stdout
(1126,16930,40)