fork(26) download 
  1. {-# OPTIONS_GHC -O2 #-}
  2. {-# LANGUAGE BangPatterns #-}
  3.  
  4. -- rosettacode.org/wiki/Hamming_numbers#Haskell    
  5.  
  6. module Main where         -- cf. ideone.com/k8PU3, GXh4P0, 01dpQu
  7.  
  8. import Data.List 
  9. import Data.Function
  10.  
  11. hamming = 1 : map (2*) hamming `union` map (3*) hamming `union` map (5*) hamming
  12.   where
  13.     union a@(x:xs) b@(y:ys) = case compare x y of
  14.             LT -> x : union  xs  b
  15.             EQ -> x : union  xs  ys
  16.             GT -> y : union  a   ys
  17.  
  18. main =
  19.  do 
  20.   s <- getLine
  21.   case s of
  22.    "a" -> do print $ take 20 hamming
  23.              print  (hamming !! 1690, hamming !! 1691)
  24.              print $ hamming !! (1000000-1)          -- 9 MB: releases the prefix of the list
  25.    "b" -> do
  26.              print $ hamming !! (1000000-1)
  27.              print $ hamming !!  1000000      -- 77 MB: does NOT release the prefix (is needed twice)
  28.    "c" -> do
  29.              mapM_ print $ take 2 $ drop 999999 hamming  -- 9 MB: used once, prefix gc-d
  30.    "d" -> do
  31.              let (_ ,(r,t))   = nthHam (1000000)         -- 4 MB: stores upper band only
  32.                  (_ ,(r2,t2)) = nthHam (1000001)
  33.              print (t,trival t)         
  34.              print (t2,trival t2)         
  35.    _   -> do                                          --    10^8:  4 MB          0.27 sec
  36.              let (nb,(r,t)) = nthHam (read s)         --    10^9:  6 MB (-4=2)   1.42 sec   ~n^0.7
  37.              print (nb, t, showLogVal 2 r)           --   10^10: 14 MB (-4=10)  7.20 sec   ~n^0.7
  38.              -- print (case t of (i,j,k)-> 2^i*3^j*5^k)
  39.  
  40. showLogVal :: Double -> Double -> String
  41. showLogVal base logval = show (10**y) ++ "E+" ++ show x
  42.    where (x,y) = properFraction (logBase 10 base * logval) 
  43.  
  44. trival (i,j,k)    = 2^i * 3^j * 5^k
  45.  
  46. -- directly find n-th Hamming number (base 1, from 2), in ~ O(n^{2/3}) time
  47. -- based on "top band" idea by Louis Klauder from DDJ discussion,
  48. -- by Will Ness, original post: drdobbs.com/blogs/architecture-and-design/228700538
  49. {-
  50. ln2 = log 2;  ln3 = log 3;  ln5 = log 5
  51. logval (i,j,k)    = fromIntegral i*ln2 + fromIntegral j*ln3 + fromIntegral k*ln5
  52. estval n          = (6*ln2*ln3*ln5* fromIntegral n)**(1/3)  -- estimated logval
  53. rngval n       
  54.     | n > 500000  = (1.698 , 0.0050)                          -- empirical
  55.     | n > 50000   = (1.693 , 0.0100)                          --  estimation
  56.     | n > 500     = (1.66  , 0.0500)                          --   correction  
  57.     | n > 1       = (1.56  , 0.2000)                          --    (dist,width)
  58.     | otherwise   = (1.56  , 0.4000) 
  59.  
  60.  -- estval(n) = log (M[n]) = ln2 * logBase 2 (M[n])
  61. -}
  62. ww = logBase 2 30 / 2
  63. -- lb3 :: Double
  64. lb3 = logBase 2 3;  lb5 = logBase 2 5   -- lb3 == log 3/log 2,  lb5 == log 5/log 2
  65. -- logval2 :: (Int,Int,Int) -> Double
  66. logval2 (i,j,k)    = fromIntegral i + fromIntegral j*lb3 + fromIntegral k*lb5
  67. -- estval2 :: Integer -> Double
  68. estval2 n          = (6*lb3*lb5* fromIntegral n)**(1/3)  -- estimated logval **Base 2**
  69. rngval2 n           --  -1/v          +2/v      seems to works too
  70.     | n > 500000  = (ww-(3/estval2 n), 6/estval2 n)      -- space tweak !!! (thx, GBG!)
  71.     | n > 500000  = (2.4496 , 0.0076 )                          -- empirical
  72.     | n > 50000   = (2.4424 , 0.0146 )                          --  estimation
  73.     | n > 500     = (2.3948 , 0.0723 )                          --   correction - base 2
  74.     | n > 1       = (2.2506 , 0.2887 )                          --    (dist,width)
  75.     | otherwise   = (2.2506 , 0.5771 )  
  76.  
  77. -- nthHam :: Integer -> ( (Int,Int), (Double, (Int, Int, Int)))  -- ( 64-bits: use Int!!    NB! )
  78. --                    _band size stuff_  -- (*1*)
  79.  
  80. nthHam :: Integer -> ( (Bool, Double, Double), (Double, (Int, Int, Int)))
  81. --                       _max logval, min delta in band_
  82. nthHam n                                                      -- n: 1-based 1,2,3...
  83.   | w >= 1   = error $ "Breach of contract: (w < 1):  " ++ show (w)
  84.   | m <  0   = error $ "Not enough triples generated: " ++ show (c,n)
  85.   | m >= nb  = error $ "Generated band is too narrow: " ++ show (m,nb)
  86.   | True     = -- ((m,nb),res) -- m=target index in band from top, nb=band's length  
  87.                ( ( isTrulySorted ,
  88.                    fst (head s),         -- (*1*)
  89.                       minimum -- By (compare `on` abs)
  90.                          (zipWith (-) (map fst s) (tail (map fst s))) )
  91.                  , res)
  92.  where 
  93.   isTrulySorted = and [a>b | let z=map (trival.snd) s, (a,b)<-zip z (tail z)]
  94.   (d,w)   = rngval2 n                                         -- correction dist, width
  95.   hi      = estval2 n - d                                     --   hi > logval2 > hi-w 
  96.   (c,b)   = foldl_                                            -- total count, the band
  97.              (\(!c,!b) (i,t) -> case t of [] -> (i+c,b)  
  98.                                           [x] -> (i+c,x:b)) 
  99.              (0::Integer,[])                                  -- ( 64bit: use Int!!!     NB! ) 
  100.              -- (sum *** concat) . unzip $
  101.              [ ( fromIntegral i+1,                            -- total triples w/ this (j,k)
  102.                  [ (r,(i,j,k)) | frac < w ] )                 -- store it, if inside band
  103.                | k <- [ 0 .. floor ( hi   /lb5) ],  let p = fromIntegral k*lb5,    
  104.                  j <- [ 0 .. floor ((hi-p)/lb3) ],  let q = fromIntegral j*lb3 + p,
  105.                  let (i,frac) = pr  (hi-q) ;            r = hi - frac  -- r = i + q 
  106.                ] where     pr = properFraction                -- pr 1.24 => (1,0.24)
  107.   (m,nb)  = ( fromIntegral $ c - n, length b )                -- m 0-based from top, |band|
  108.   (s,res) = ( sortBy (flip compare `on` fst) b, s!!m )        -- sorted decreasing, result
  109.   foldl_  = foldl'
Success #stdin #stdout 3.26s 6924KB
comments (?)
stdin 
100000000000

c  0.85s-8.0MB
a  0.87s-8.0MB 

is Double precision enough?
1T:   (28052.292341476037,2.983142621815205e-9) 14 significant digits needed.   iffy... but OK
100B: (13019.406212824639,4.964022082276642e-9) 14 significant digits needed. 
10B:   (6041.758782302383,2.292654244229198e-8) 12 significant digits needed. 
1B:    (2803.022191612377,1.325179255218245e-7) 11 significant digits needed.  


1T   (1126,16930, 40), "3.814325005007765E+8444"            ( GXh4P0 20% faster,  
100B (10178,1384,279), "1.705075482875041E+3919"                with explicit loops)
10B  (2177,  8, 1659), "5.629992181012814E+1818" 
1B   (1334, 335, 404), "6.216075755565590E+843"  
100M (  2,  454, 249), "1.814014330961032E+391"
____
TIMINGS, without the isTrulySorted calculation which slowes it down by ~20% for 100B:
----
             target idx    
            (  from top,band size ) 
1T   13.39s-7.9MB (8609,22872)   37.63%    ~ n^0.69 time, n^0.33 space
100B  2.74s-6.9MB (3993,10609)   37.64%          --------    (no tweak: ~3.2s; ~175,000 band size)
50B   1.66s-6.9MB (3167,8421 )   37.61%             |
10B   0.49s-5.9MB (1854,4928 )   37.62%             n^(1/3) -- band size
1B    0.05s-4.9MB ( 862,2288 )   37.67%             |
100M  0.01s-4.9MB ( 399,1061 )   37.61%          --------
10M   0.00s-4.9MB ( 183,492  )   37.20% 
5M    0.00s-4.9MB ( 150,391  )   38.36%     ( no tweak:  5M 0-4.9  (74,238) 0.311%  )
1M    0.00s-4.9MB (  86,231  )   37.23%     ( no tweak:  1M 0-4.9  (12,79)  0.1519% )
stdout 
((True,13019.406212824639,4.964022082276642e-9),(10178,1384,279),"1.7050754828750412E+3919")
https://ideone.com/q3fma
http://stackoverflow.com/a/10160054/849891
language:
Haskell (ghc 8.4.4)
created:
8 years ago
visibility:
public

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!