fork download
  1. {-# OPTIONS -O3 -XStrict -XBangPatterns #-}
  2.  
  3. import Data.List (sortBy, foldl')
  4. import Data.Function (on)
  5.  
  6. main = let (r,t) = nthHam 1000000000000 in print t -- >> print (trival t)
  7.  
  8. trival (i,j,k) = 2^i * 3^j * 5^k
  9.  
  10. nthHam :: Int -> (Double, (Int, Int, Int)) -- ( 64bit: use Int!!! NB! )
  11. nthHam n -- n: 1-based: 1,2,3...
  12. | n <= 0 = error $ "n is 1--based: must be n > 0: " ++ show n
  13. | n < 2 = ( 0.0, (0, 0, 0) ) -- trivial case so estimation works for rest
  14. | w >= 1 = error $ "Breach of contract: (w < 1): " ++ show w
  15. | m < 0 = error $ "Not enough triples generated: " ++ show ((c,n) :: (Int, Int))
  16. | m >= nb = error $ "Generated band is too narrow: " ++ show (m,nb)
  17. | otherwise = sortBy (flip compare `on` fst) b !! m -- m-th from top in sorted band
  18. where
  19. lb3 = logBase 2 3; lb5 = logBase 2 5; lb30_2 = logBase 2 30 / 2
  20. v = (6*lb3*lb5* fromIntegral n)**(1/3) - lb30_2 -- estimated logval, base 2
  21. estval n = (v + (1/v), 2/v) -- the space tweak! (thx, GBG!)
  22. (hi,w) = estval n -- hi > logval > hi-w
  23. m = fromIntegral (c - n) -- target index, from top
  24. nb = length (b :: [(Double, (Int, Int, Int))]) -- length of the band
  25. (c,b) = foldl_ (\(c,b) (i,t)-> let c2=c+i in c2 `seq` -- ( total count, the band )
  26. case t of []-> (c2,b);[v]->(c2,v:b) ) (0,[]) -- ( =~= mconcat )
  27. [ ( fromIntegral i+1, -- total triples w/ this (j,k)
  28. [ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band
  29. | k <- [ 0 .. floor ( hi /lb5) ], let p = fromIntegral k*lb5,
  30. j <- [ 0 .. floor ((hi-p)/lb3) ], let q = fromIntegral j*lb3 + p,
  31. let (i,frac) = pr (hi-q) ; r = hi - frac -- r = i + q
  32. ] where pr = properFraction -- pr 1.24 => (1,0.24)
  33. foldl_ = foldl'
  34.  
Success #stdin #stdout 0.78s 0KB
stdin
Standard input is empty
stdout
(1126,16930,40)