{-# OPTIONS_GHC -O2 #-} {-# LANGUAGE BangPatterns #-} -- rosettacode.org/wiki/Hamming_numbers#Haskell module Main where -- cf. ideone.com/k8PU3, GXh4P0, 01dpQu import Data.List import Data.Function hamming = 1 : map (2*) hamming `union` map (3*) hamming `union` map (5*) hamming where union a@(x:xs) b@(y:ys) = case compare x y of LT -> x : union xs b EQ -> x : union xs ys GT -> y : union a ys main = do s <- getLine case s of "a" -> do print $ take 20 hamming print (hamming !! 1690, hamming !! 1691) print $ hamming !! (1000000-1) -- 9 MB: releases the prefix of the list "b" -> do print $ hamming !! (1000000-1) print $ hamming !! 1000000 -- 77 MB: does NOT release the prefix (is needed twice) "c" -> do mapM_ print $ take 2 $ drop 999999 hamming -- 9 MB: used once, prefix gc-d "d" -> do let (_ ,(r,t)) = nthHam (1000000) -- 4 MB: stores upper band only (_ ,(r2,t2)) = nthHam (1000001) print (t,trival t) print (t2,trival t2) _ -> do -- 10^8: 4 MB 0.27 sec let (nb,(r,t)) = nthHam (read s) -- 10^9: 6 MB (-4=2) 1.42 sec ~n^0.7 print (nb, t, showLogVal 2 r) -- 10^10: 14 MB (-4=10) 7.20 sec ~n^0.7 -- print (case t of (i,j,k)-> 2^i*3^j*5^k) showLogVal :: Double -> Double -> String showLogVal base logval = show (10**y) ++ "E+" ++ show x where (x,y) = properFraction (logBase 10 base * logval) trival (i,j,k) = 2^i * 3^j * 5^k -- directly find n-th Hamming number (base 1, from 2), in ~ O(n^{2/3}) time -- based on "top band" idea by Louis Klauder from DDJ discussion, -- by Will Ness, original post: drdobbs.com/blogs/architecture-and-design/228700538 {- ln2 = log 2; ln3 = log 3; ln5 = log 5 logval (i,j,k) = fromIntegral i*ln2 + fromIntegral j*ln3 + fromIntegral k*ln5 estval n = (6*ln2*ln3*ln5* fromIntegral n)**(1/3) -- estimated logval rngval n | n > 500000 = (1.698 , 0.0050) -- empirical | n > 50000 = (1.693 , 0.0100) -- estimation | n > 500 = (1.66 , 0.0500) -- correction | n > 1 = (1.56 , 0.2000) -- (dist,width) | otherwise = (1.56 , 0.4000) -- estval(n) = log (M[n]) = ln2 * logBase 2 (M[n]) -} ww = logBase 2 30 / 2 -- lb3 :: Double lb3 = logBase 2 3; lb5 = logBase 2 5 -- lb3 == log 3/log 2, lb5 == log 5/log 2 -- logval2 :: (Int,Int,Int) -> Double logval2 (i,j,k) = fromIntegral i + fromIntegral j*lb3 + fromIntegral k*lb5 -- estval2 :: Integer -> Double estval2 n = (6*lb3*lb5* fromIntegral n)**(1/3) -- estimated logval **Base 2** rngval2 n -- -1/v +2/v seems to works too | n > 500000 = (ww-(3/estval2 n), 6/estval2 n) -- space tweak !!! (thx, GBG!) | n > 500000 = (2.4496 , 0.0076 ) -- empirical | n > 50000 = (2.4424 , 0.0146 ) -- estimation | n > 500 = (2.3948 , 0.0723 ) -- correction - base 2 | n > 1 = (2.2506 , 0.2887 ) -- (dist,width) | otherwise = (2.2506 , 0.5771 ) -- nthHam :: Integer -> ( (Int,Int), (Double, (Int, Int, Int))) -- ( 64-bits: use Int!! NB! ) -- _band size stuff_ -- (*1*) nthHam :: Integer -> ( (Bool, Double, Double), (Double, (Int, Int, Int))) -- _max logval, min delta in band_ nthHam n -- n: 1-based 1,2,3... | w >= 1 = error $ "Breach of contract: (w < 1): " ++ show (w) | m < 0 = error $ "Not enough triples generated: " ++ show (c,n) | m >= nb = error $ "Generated band is too narrow: " ++ show (m,nb) | True = -- ((m,nb),res) -- m=target index in band from top, nb=band's length ( ( isTrulySorted , fst (head s), -- (*1*) minimum -- By (compare `on` abs) (zipWith (-) (map fst s) (tail (map fst s))) ) , res) where isTrulySorted = and [a>b | let z=map (trival.snd) s, (a,b)<-zip z (tail z)] (d,w) = rngval2 n -- correction dist, width hi = estval2 n - d -- hi > logval2 > hi-w (c,b) = foldl_ -- total count, the band (\(!c,!b) (i,t) -> case t of [] -> (i+c,b) [x] -> (i+c,x:b)) (0::Integer,[]) -- ( 64bit: use Int!!! NB! ) -- (sum *** concat) . unzip $ [ ( fromIntegral i+1, -- total triples w/ this (j,k) [ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band | k <- [ 0 .. floor ( hi /lb5) ], let p = fromIntegral k*lb5, j <- [ 0 .. floor ((hi-p)/lb3) ], let q = fromIntegral j*lb3 + p, let (i,frac) = pr (hi-q) ; r = hi - frac -- r = i + q ] where pr = properFraction -- pr 1.24 => (1,0.24) (m,nb) = ( fromIntegral $ c - n, length b ) -- m 0-based from top, |band| (s,res) = ( sortBy (flip compare `on` fst) b, s!!m ) -- sorted decreasing, result foldl_ = foldl'