fork(4) download
copy 
  1. {-# OPTIONS_GHC -O2 #-}
  2. {-# LANGUAGE BangPatterns #-}
  3.  
  4. -- http://r...content-available-to-author-only...e.org/wiki/Hamming_numbers#Haskell
  5.  
  6. module Main where
  7. import Data.List 
  8. import Data.Function
  9.  
  10. hamming = 1 : map (2*) hamming `union` map (3*) hamming `union` map (5*) hamming
  11.   where
  12.     union a@(x:xs) b@(y:ys) = case compare x y of
  13.             LT -> x : union  xs  b
  14.             EQ -> x : union  xs  ys
  15.             GT -> y : union  a   ys
  16.  
  17. main = do 
  18.   s <- getLine
  19.   case s of
  20.    "1" -> do
  21.              print $ take 20 hamming
  22.              print  (hamming !! 1690, hamming !! 1691)
  23.              print $ hamming !! (1000000-1)          -- 9 MB: releases the prefix of the list
  24.    "2" -> do
  25.              print $ hamming !! (1000000-1)
  26.              print $ hamming !!  1000000             -- 77 MB: does NOT release the prefix (is needed twice)
  27.    "11" -> do
  28.              let (r,t) = nthHam (1000000-1)          -- 4 MB: stores upper band only
  29.              print (t,trival t)         
  30.    _   -> do        
  31.              let (r,t) = nthHam (read s)             --    10^8: 5 MB   0.29 sec
  32.              print t                                 --    10^9: 6 MB   1.52 sec   O(n^0.7)
  33.  
  34. -- directly find n-th Hamming number (base 1, from 2), in ~ O(n^{2/3}) time
  35. -- by Will Ness, based on "band" idea by Louis Klauder from DDJ discussion
  36. --   http://d...content-available-to-author-only...s.com/blogs/architecture-and-design/228700538
  37.  
  38. ln2 = log 2;  ln3 = log 3;  ln5 = log 5
  39. logval (i,j,k)    = fromIntegral i*ln2 + fromIntegral j*ln3 + fromIntegral k*ln5
  40. trival (i,j,k)    = 2^i * 3^j * 5^k
  41. estval n          = (6*ln2*ln3*ln5* fromIntegral n)**(1/3)  -- estimated logval
  42. rngval n       
  43.     | n > 500000  = (1.698 , 0.0050)                          -- empirical
  44.     | n > 50000   = (1.693 , 0.0100)                          --  estimation
  45.     | n > 500     = (1.66  , 0.0500)                          --   correction  
  46.     | n > 1       = (1.56  , 0.2000)                          --    (dist,width)
  47.     | otherwise  = (1.56  , 0.4000)  
  48.  
  49. nthHam n1                                                     -- n1: 1-based 2,3...
  50.   | w >= ln2 = error $ "Breach of contract: w >= ln2  " ++ show (w,ln2)
  51.   | m <  0   = error $ "Not enough triples generated: " ++ show (c,n)
  52.   | m >= nb  = error $ "Generated band is too narrow: " ++ show (m,nb)
  53.   | True     = res 
  54.  where 
  55.   n       = n1 + 1                                            -- n: 1-based 1,2,3...
  56.   (d,w)   = rngval n                                          -- correction dist, width
  57.   (w2,hi) = ( w/ln2, estval n - d )                           --   hi > logval > hi-w 
  58.   (c,b)   = foldl'                                            -- total count, the band
  59.              (\(!c,!b) (i,t) -> case t of [] -> (i+c,b)
  60.                                           [x] -> (i+c,x:b)) 
  61.              (0,[])
  62.              [ ( i+1,                                         -- total triples w/ this (j,k)
  63.                  [ (r,(i,j,k)) | frac < w2 ] )                -- store it, if inside band
  64.                | k <- [ 0 .. floor ( hi   /ln5) ],  let p = fromIntegral k*ln5,    
  65.                  j <- [ 0 .. floor ((hi-p)/ln3) ],  let q = fromIntegral j*ln3 + p,
  66.                  let (i,frac) = pr ((hi-q)/ln2) ;        r = fromIntegral i*ln2 + q 
  67.                ] where     pr = properFraction
  68.   (m,nb)  = ( fromInteger $ c - n, length b )               -- m 0-based from top, |band|
  69.   (s,res) = ( sortBy (flip compare `on` fst) b, s!!m )      -- sorted decreasing, result
stdin
copy 
1
compilation info 
[1 of 1] Compiling Main             ( prog.hs, prog.o )
Linking prog ...
stdout
copy 
[1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36]
(2125764000,2147483648)
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
https://ideone.com/PKxza
language:
Haskell (ghc 8.4.4)
created:
14 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!