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