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

to print the full value as well, uncomment one of the
  `-- print (case t of (i,j,k) -> 2^i*3^j*5^k)`
lines in `main`, at its very end

------
2020-03-24: use `Int` in nthHam, now on 64 bit, for speed:

1B  0.02s 6.0MB
((True,2803.022191612378,1.325179255218245e-7),(1334,335,404),"6.216075755562335E+843")
1T  2.01s 15.1MB
((True,28052.292341476037,2.979504643008113e-9),(1126,16930,40),"3.814325005007765E+8444")
4T  5.81s 22.2MB   16 digits used.... 
                  still good (as evidenced by the `True` below), but really pushing it.
((True,44531.67942695671,7.275957614183426e-11),(16348,16503,873),"2.3509832704071347E+13405")
10T   11.13s 26.4MB
((True,60439.66391703062,7.275957614183426e-11),(18187,23771,1971),"1.4182592613214273E+18194")
13T   14.44s 30.4MB    ...still good
((True,65963.64327239885,5.820766091346741e-11),(28648,21308,1526),"1.0845209980457384E+19857")

same code on tio:
10T   16.77s
35T   38.84s 
((True,91766.48002166452,5.820766091346741e-11),(13824,2133,32112),"2.904528312044042E+27624")
60T   53.76s
((True,109828.1626252908,5.820766091346741e-11),(28123,45913,3848),"3.7265992051315013E+33061")
65T   53.98s
((True,112797.98518694332,5.820766091346741e-11),(26349,21177,22776),"3.7755976376763294E+33955")
68T   55.88s
((True,114507.3424180621,5.820766091346741e-11),(28179,16558,25877),"1.395679394278755E+34470")
70T   59.57s
((True,115619.15754088841,5.820766091346741e-11),(13125,13687,34799),"6.831047840502602E+34804")

on home machine:
100T: 368.13s     n^0.708
((True,130216.14085518729,5.820766091346741e-11),(88324,876,17444),"9.211106366443564E+39198")

140T: 466.69s     n^0.705
((True,145671.64804685948,5.820766091346741e-11),(9918,24002,42082),"3.432221006522245E+43851")

170T: 383.26s     n^0.076   ---FAULTY---
((False,155411.25012054382,0.0),(77201,27980,14584),"2.8050819114774215E+46783")
------

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"
____
old (32-bit) 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
copy 
((True,28052.292341476037,2.979504643008113e-9),(1126,16930,40),"3.814325005007765E+8444")
381432500501368203378087989364684566003196981728287861209871741226226525212721355737268953258289365867404981597404054421789109760844500952247925146099494880428508529216246244511244867744925683802530879101066553238177517252767636184594334387309922212648235124288475749806206347656448649307449094910179817154595317474803832693242133737544677246725893557633654926313190236770066360622938574559713516391564144547396432295444647508095324602670011327673727739957987675308549713969242190198751614368654552970058346461965230959479549183120034812416513902203057203906014674293111023568658682172646478061761547592733863996763247542352604485351550500259020895056662371226777511522717096663981648958925160590781203849840250387080913587872019962833628471286020944446036116227561664422800434640303596225654342128319942334560763248219521958112980122656154322288782197044471919701456432396711460047592839239991995065396986351375936676582038764364135956974830358423813236898274231178796818725438341925717416647797193232651544867460829679476323177007770129492717496915465821550889348573821649618385759888429506048769746984318779379387235508725017626931291783495454078633856251352965367975156315500138512273903936359696710772951672000764780797534135617532615069653722327470053996537987966507375483235828950483009204942019944815486836876541835728053944013919256953575734463016041957802358498736651342458150091298395471511454596373412462056385448933043835369209570333091893485364162359777824091047663292556531848442316429641804084227404418581764178118799902088490330934415783419436766111289751728080616887467259380495593993999007416205941205205648099105889610584749882260704486606343851508816670299483639904252062724303543588169194142282157474210801912915786508784824094915984977696488857697688837471802181949631171427892810191259844779496613305824486287004098713333319920755014545092234973161765100732819317337968563808224916203530649482072304780251606930250679585076842569969013704981775998788013146577781010855892536687825783816873026869460088122790746302240414681660795730017210102260941819954392050654877165897092829136349360646270693828793890131656110085555102549404684366912948435560177729319213268641684600573665767060907665495486646622569960715616755051756190976734484896258475838165360541628697955550814600563240253989152549004757785712811084127267978889840332750151081858425794804819882698588743364095268350947868326169540540927401161514106448366645429921763405608340439515742058526616878358948974369782609855655675294824544796287426791518372131265927196232745663598948967696042989478611802894904236740139363217489156764375914289396097600245602455102665784834589391013700595861674262694959296141311362268034856325830685732072966898388781459190722295245407456302472041305892463709388574440918398498912456769938061673445888390451220594643942529155114279393017006570423565313506923730210404246978328087096083105528955910796927360711527123422727091919503910519570848679688959076293500742647850023183903147481723829114314447574070129805594179742737206447355799817808476580365839278299208478092312652752038210700664065535500101795705182938559765635161132444138452049446242866671947399765832117854684941405518872019649121058590993958524773370229558868931306124971557064586217966030355901338403646166993476691398603434729560312281007995073275898516459971220715657699387523541330373102873040610222293381440297359855604813713487045065173446504818338471581488928367463848756107869477386616119827264615280682634297173123261177373540526229015339553181728486997297652385088193337393105864888294893549092590589234315492057527958083663610553788816353890227097781381357352939605796261046157741589539421216828857527824381073392141187857645455838419363948924419249789707529869828445816097149777574457846168863021576195719281270085256789215456845232276418086131573167048399735999446047308893077404474069962516873750976531134605177381860006560751623956248369842277829338354338239593657154667339823706858657828846076357564141041144132653918540614567364428999424953532960255704726380063206414477874473011744905086193323212196230461016234776009019556145289115469651975034868214875465621658690986611676250782377764355946861847552708036446068853697203543193591291479535305020014678268365390166542449626331277902891960154198320640447275407283842209433526300227914464028822306278110016143575327647214340233047351066244266576781395298905370583147237169346233738497514110562281408609355349560048891725146122806891273216369639851077436406592893753729690835960446164271693935812602839579229273571311436712152831655030039235186756793939086768598572119711208661869571225003654933600358235541115015140765308202890492006769502831942411647119268717661579552940663489637401616797752823585038415670347522491030842381108052530483847642586054403601447134420408642823870164041381659369018513999163584560593071103729523430358336622667039689019463450035988415451651556138276294562084712233906409886132878792332009910602861356452290649901406349299580539448017304889428440240802135729922126748998968758135863691971362531429661032283162834322098678183994246815437051118481713336003797985203440751253618910774614857932245930438370672346152009783488939557110237315572630011188807032434684323410978888350161084958463242054583107168129809365906510277520834199584853707564975390831892479367377747704598145178429542726047766484073513798982629254881015589607774241895593929484591210459733972089725801226712669526167643559949974875548078686724927959894386937328668387305412202357110352655277305103854151248777748003912914187166760456241514381336238622922176322074649573786643542323409287130562729156125057560479177108923976620833793789796446288379720457136250236958013583837301292182469558879831232003762891433765287723404287606949817910100748412691304937225478604543944712581388884826063975760055469108001222120792097236141905327681074374503535383573887538394160004316349174529760474207668106657055423044749824787960307145882371169230102374074351105763898972107677525315719770924516050502350592454339732826134162740274248974361445949386533845057455877594711936404651222213700068778072085053124279558647251710675496148419947187710144483465061491027072800457529151879293653843768514232911643040932812315558121776764357308509623935085751986865128160332811976103135474624261041393051559430318488475203773394340924758715006468013654729377120628769731870289123185162241978318947284851186423539183667271806058237944470875029035834532404348499800745036015383845562794552147701458176737502147741344932704133451193644858265528208740292334938649754442528229593081859268043589651263172114838672296433455437884811594947448395870528675891271525327850454305478319284798674414283610758455864692025010128412213817083204073411513858208943478155840652861722720103650369636812940260428918875534584092687774013464914555586538728185037492887768408661425559897617085879377999322482768408377153456578353192600088691156247278285376650302918020364100406696826182661984739230472483569772685421656899496425201205051292520579762192167734835433528278166901986761398474772205192235112447073768501872380291313256166658437838705840725910673067001694599268891466041051301357758142028926221241150737563383094473077416353682405442288466206749465489838345288022292312683862458722942463812695992417502575991713081094596966691251333855046350205503351216355174898522776275164001673042505352494099356428324565551222456829595241110281915073149143722112106597387867484881077143415016534201519380037945223589464410350907465343341743447452315897697439240008325903230535647680590852852742770572098939348138339728240185339597903790459538131474918499054223109084798550228061965714109171636789980604603951414783253430837758215397191494180722955449120576644783948206405282461782774370469388077935264165516109137201074028611700340924754761566309315919658543055145034524177191198831130692839923393256394546847019892233823877166696730848415975077044614276994220212117518343302008536715617994671386526822378431167884328996507773315531097327141735084258815563626491801736919645485414765901463436227902736185065138747561652147253349857191230736498725073474761761827043174919382521944349056386305186481198173455293985211104739115063577688592953395294676810655962894343645229199990527683439765973171737546615533185149900878528865499358951458942761341111026253722935261672690109408935429704509667276178562289568318381821903503360000000000000000000000000000000000000000
https://ideone.com/gHPmTo
http://stackoverflow.com/a/10160054/849891
language:
Haskell (ghc 8.4.4)
created:
4 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!