{-# OPTIONS_GHC -O2 #-}
{-# LANGUAGE BangPatterns #-}
-- rosettacode.org/wiki/Hamming_numbers#Haskell -- derived from ideone.com/q3fma
module Main where
import Data.List
import Data.Function
newtype Ham
= Ham
(Double,(Int,Int,Int)) -- 1M: OK: ("5.193127804483973E+83",(55,47,64))
Ham (_,a) == Ham (_,b) = a == b
show (Ham
(x
,(i
,j
,k
))) = show (showLogVal
2 x
,(i
,j
,k
)) -- (x,2^i*3^j*5^k)
showLogVal base logval
= show (10**y
) ++ "E+" ++ show x
mul2 (Ham (_,(b,c,d))) = Ham ((_I b+1.0) + _I c *lb3 + _I d *lb5, (b+1,c,d))
mul3 (Ham (_,(b,c,d))) = Ham ( _I b + (_I c+1.0)*lb3 + _I d *lb5, (b,c+1,d))
mul5 (Ham (_,(b,c,d))) = Ham ( _I b + _I c *lb3 + (_I d+1.0)*lb5, (b,c,d+1))
one = Ham (0.0,(0,0,0))
hamming
= one :
map mul2 hamming `union`
map mul3 hamming `union`
map mul5 hamming
-- [Ham] -- 1 : map (2*) hamming `union` map (3*) hamming `union` map (5*) hamming -- [Integer]
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
case s of
('x':s
) -> do print $ hamming
!! (read s
-1) -- 1-based: 1-st is (0.0,(0,0,0)) print (hamming
!! 1690, hamming
!! 1691) print $ hamming
!! (1000000-1) -- 9 MB: releases the prefix of the list "b" -> do
print $ hamming
!! (1000000-1) -- 77 MB: does NOT release the prefix print $ hamming
!! 1000000 -- (is needed twice) "c" -> do
"d" -> do
let (r,t) = nthHam (1000000) -- 4 MB: stores upper band only
(r2,t2) = nthHam (1000001)
_ -> do -- 10^8: 4 MB 0.27 sec
let (r
,t
) = nthHam
(read s
) -- 10^9: 6 MB (-4=2) 1.42 sec O(n^0.7) print t
-- 10^10: 14 MB (-4=10) 7.20 sec O(n^0.7)
-- 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
trival (i,j,k) = 2^i * 3^j * 5^k
estval n
= (6*ln2
*ln3
*ln5
* fromIntegral n
)**(1/3) -- estimated logvalrngval 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)
-- estval(n) = log (M[n]) = ln2 * logBase 2 (M[n])
lb3
= logBase 2 3; lb5
= logBase 2 5 -- lb3 == log 3/log 2, lb5 == log 5/log 2estval2 n
= (6*lb3
*lb5
* fromIntegral n
)**(1/3) -- estimated logval **Base 2**rngval2 n
| 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)
nthHam n -- n: 1-based 1,2,3...
| w
>= 1 = error $ "Breach of contract (w < 1): " ++ show (w
) | m
< 0 = error $ "The estimate is too small: " ++ show (c
,n
) | m
>= nb
= error $ "Generated band is too narrow: " ++ show (m
,nb
) | True = res
where
(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,[])
[ ( i+1, -- total triples w/ this (j,k)
[ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band
let (i,frac) = pr (hi-q) ; r = hi - frac -- r = i + q
(s
,res
) = ( sortBy
(flip compare `on`
fst) b
, s
!!m
) -- sorted decreasing, result
ey0jIE9QVElPTlNfR0hDIC1PMiAjLX0Key0jIExBTkdVQUdFIEJhbmdQYXR0ZXJucyAjLX0KCi0tIHJvc2V0dGFjb2RlLm9yZy93aWtpL0hhbW1pbmdfbnVtYmVycyNIYXNrZWxsICAgICAgLS0gZGVyaXZlZCBmcm9tIGlkZW9uZS5jb20vcTNmbWEKCm1vZHVsZSBNYWluIHdoZXJlCmltcG9ydCBEYXRhLkxpc3QgCmltcG9ydCBEYXRhLkZ1bmN0aW9uCgpuZXd0eXBlIEhhbSA9IEhhbSAoRG91YmxlLChJbnQsSW50LEludCkpICAtLSAxTTogT0s6ICgiNS4xOTMxMjc4MDQ0ODM5NzNFKzgzIiwoNTUsNDcsNjQpKQoKaW5zdGFuY2UgRXEgSGFtIHdoZXJlCiAgSGFtIChfLGEpID09IEhhbSAoXyxiKSA9IGEgPT0gYgoKaW5zdGFuY2UgT3JkIEhhbSB3aGVyZQogIGNvbXBhcmUgKEhhbSAoYSxfKSkgKEhhbSAoYixfKSkgPSBjb21wYXJlIGEgYgoKaW5zdGFuY2UgU2hvdyBIYW0gd2hlcmUKICBzaG93IChIYW0oeCwoaSxqLGspKSkgPSBzaG93IChzaG93TG9nVmFsIDIgeCwoaSxqLGspKSAgIC0tICh4LDJeaSozXmoqNV5rKQogIApzaG93TG9nVmFsIGJhc2UgbG9ndmFsID0gc2hvdyAoMTAqKnkpICsrICJFKyIgKysgc2hvdyB4CiAgIHdoZXJlICh4LHkpID0gcHJvcGVyRnJhY3Rpb24gKGxvZ0Jhc2UgMTAgYmFzZSAqIGxvZ3ZhbCkgCgpfSSBpID0gZnJvbUludGVncmFsIGkKCm11bDIgKEhhbSAoXywoYixjLGQpKSkgPSBIYW0gKChfSSBiKzEuMCkgKyAgX0kgYyAgICAgKmxiMyArICBfSSBkICAgICAqbGI1LCAoYisxLGMsZCkpCm11bDMgKEhhbSAoXywoYixjLGQpKSkgPSBIYW0gKCBfSSBiICAgICAgKyAoX0kgYysxLjApKmxiMyArICBfSSBkICAgICAqbGI1LCAoYixjKzEsZCkpCm11bDUgKEhhbSAoXywoYixjLGQpKSkgPSBIYW0gKCBfSSBiICAgICAgKyAgX0kgYyAgICAgKmxiMyArIChfSSBkKzEuMCkqbGI1LCAoYixjLGQrMSkpCgpvbmUgPSBIYW0gKDAuMCwoMCwwLDApKQoKaGFtbWluZyA9IG9uZSA6IG1hcCBtdWwyIGhhbW1pbmcgYHVuaW9uYCBtYXAgbXVsMyBoYW1taW5nIGB1bmlvbmAgbWFwIG11bDUgaGFtbWluZyAgIC0tIFtIYW1dCiAgICAgICAtLSAxIDogbWFwICgyKikgaGFtbWluZyBgdW5pb25gIG1hcCAoMyopIGhhbW1pbmcgYHVuaW9uYCBtYXAgKDUqKSBoYW1taW5nICAgICAgLS0gW0ludGVnZXJdCiAgd2hlcmUKICAgIHVuaW9uIGFAKHg6eHMpIGJAKHk6eXMpID0gY2FzZSBjb21wYXJlIHggeSBvZgogICAgICAgICAgICBMVCAtPiB4IDogdW5pb24gIHhzICBiCiAgICAgICAgICAgIEVRIC0+IHggOiB1bmlvbiAgeHMgIHlzCiAgICAgICAgICAgIEdUIC0+IHkgOiB1bmlvbiAgYSAgIHlzCgptYWluID0gZG8gCiAgcyA8LSBnZXRMaW5lCiAgY2FzZSBzIG9mCiAgICgneCc6cykgLT4gZG8gcHJpbnQgJCBoYW1taW5nICEhIChyZWFkIHMtMSkgICAgICAgLS0gMS1iYXNlZDogMS1zdCBpcyAoMC4wLCgwLDAsMCkpCiAgICJhIiAtPiBkbyBwcmludCAkIHRha2UgMjAgaGFtbWluZyAKICAgICAgICAgICAgIHByaW50ICAoaGFtbWluZyAhISAxNjkwLCBoYW1taW5nICEhIDE2OTEpCiAgICAgICAgICAgICBwcmludCAkIGhhbW1pbmcgISEgKDEwMDAwMDAtMSkgICAgICAgICAgLS0gOSBNQjogcmVsZWFzZXMgdGhlIHByZWZpeCBvZiB0aGUgbGlzdAogICAiYiIgLT4gZG8KICAgICAgICAgICAgIHByaW50ICQgaGFtbWluZyAhISAoMTAwMDAwMC0xKSAgICAgICAgICAtLSA3NyBNQjogZG9lcyBOT1QgcmVsZWFzZSB0aGUgcHJlZml4CiAgICAgICAgICAgICBwcmludCAkIGhhbW1pbmcgISEgIDEwMDAwMDAgICAgICAgICAgICAgLS0gICAgICAgICAgICAgICAgICAgIChpcyBuZWVkZWQgdHdpY2UpCiAgICJjIiAtPiBkbwogICAgICAgICAgICAgbWFwTV8gcHJpbnQgJCB0YWtlIDIgJCBkcm9wIDk5OTk5OSBoYW1taW5nICAtLSA5IE1COiB1c2VkIG9uY2UsIHByZWZpeCBnYy1kCiAgICJkIiAtPiBkbwogICAgICAgICAgICAgbGV0IChyLHQpICAgPSBudGhIYW0gKDEwMDAwMDApICAgICAgICAgICAtLSA0IE1COiBzdG9yZXMgdXBwZXIgYmFuZCBvbmx5CiAgICAgICAgICAgICAgICAgKHIyLHQyKSA9IG50aEhhbSAoMTAwMDAwMSkKICAgICAgICAgICAgIHByaW50ICh0LHRyaXZhbCB0KSAgICAgICAgIAogICAgICAgICAgICAgcHJpbnQgKHQyLHRyaXZhbCB0MikgICAgICAgICAKICAgXyAgIC0+IGRvICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tICAgIDEwXjg6ICA0IE1CICAgICAgICAgIDAuMjcgc2VjCiAgICAgICAgICAgICBsZXQgKHIsdCkgPSBudGhIYW0gKHJlYWQgcykgICAgICAgICAgICAtLSAgICAxMF45OiAgNiBNQiAoLTQ9MikgICAxLjQyIHNlYyAgTyhuXjAuNykKICAgICAgICAgICAgIHByaW50IHQgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tICAgMTBeMTA6IDE0IE1CICgtND0xMCkgIDcuMjAgc2VjICBPKG5eMC43KQoKLS0gZGlyZWN0bHkgZmluZCBuLXRoIEhhbW1pbmcgbnVtYmVyIChiYXNlIDEsIGZyb20gMiksIGluIH4gTyhuXnsyLzN9KSB0aW1lCi0tIGJhc2VkIG9uICJ0b3AgYmFuZCIgaWRlYSBieSBMb3VpcyBLbGF1ZGVyIGZyb20gRERKIGRpc2N1c3Npb24sCi0tIGJ5IFdpbGwgTmVzcywgb3JpZ2luYWwgcG9zdDogZHJkb2Jicy5jb20vYmxvZ3MvYXJjaGl0ZWN0dXJlLWFuZC1kZXNpZ24vMjI4NzAwNTM4CgpsbjIgPSBsb2cgMjsgIGxuMyA9IGxvZyAzOyAgbG41ID0gbG9nIDUKbG9ndmFsIChpLGosaykgICAgPSBmcm9tSW50ZWdyYWwgaSpsbjIgKyBmcm9tSW50ZWdyYWwgaipsbjMgKyBmcm9tSW50ZWdyYWwgaypsbjUKdHJpdmFsIChpLGosaykgICAgPSAyXmkgKiAzXmogKiA1XmsKZXN0dmFsIG4gICAgICAgICAgPSAoNipsbjIqbG4zKmxuNSogZnJvbUludGVncmFsIG4pKiooMS8zKSAgLS0gZXN0aW1hdGVkIGxvZ3ZhbApybmd2YWwgbiAgICAgICAKICAgIHwgbiA+IDUwMDAwMCAgPSAoMS42OTggLCAwLjAwNTApICAgICAgICAgICAgICAgICAgICAgICAgICAtLSBlbXBpcmljYWwKICAgIHwgbiA+IDUwMDAwICAgPSAoMS42OTMgLCAwLjAxMDApICAgICAgICAgICAgICAgICAgICAgICAgICAtLSAgZXN0aW1hdGlvbgogICAgfCBuID4gNTAwICAgICA9ICgxLjY2ICAsIDAuMDUwMCkgICAgICAgICAgICAgICAgICAgICAgICAgIC0tICAgY29ycmVjdGlvbiAgCiAgICB8IG4gPiAxICAgICAgID0gKDEuNTYgICwgMC4yMDAwKSAgICAgICAgICAgICAgICAgICAgICAgICAgLS0gICAgKGRpc3Qsd2lkdGgpCiAgICB8IG90aGVyd2lzZSAgPSAoMS41NiAgLCAwLjQwMDApIAoKIC0tIGVzdHZhbChuKSA9IGxvZyAoTVtuXSkgPSBsbjIgKiBsb2dCYXNlIDIgKE1bbl0pCgpsYjMgPSBsb2dCYXNlIDIgMzsgIGxiNSA9IGxvZ0Jhc2UgMiA1ICAgLS0gbGIzID09IGxvZyAzL2xvZyAyLCAgbGI1ID09IGxvZyA1L2xvZyAyCmxvZ3ZhbDIgKGksaixrKSAgICA9IGZyb21JbnRlZ3JhbCBpICsgZnJvbUludGVncmFsIGoqbGIzICsgZnJvbUludGVncmFsIGsqbGI1CmVzdHZhbDIgbiAgICAgICAgICA9ICg2KmxiMypsYjUqIGZyb21JbnRlZ3JhbCBuKSoqKDEvMykgIC0tIGVzdGltYXRlZCBsb2d2YWwgKipCYXNlIDIqKgpybmd2YWwyIG4gICAgICAgCiAgICB8IG4gPiA1MDAwMDAgID0gKDIuNDQ5NiAsIDAuMDA3NiApICAgICAgICAgICAgICAgICAgICAgICAgICAtLSBlbXBpcmljYWwKICAgIHwgbiA+IDUwMDAwICAgPSAoMi40NDI0ICwgMC4wMTQ2ICkgICAgICAgICAgICAgICAgICAgICAgICAgIC0tICBlc3RpbWF0aW9uCiAgICB8IG4gPiA1MDAgICAgID0gKDIuMzk0OCAsIDAuMDcyMyApICAgICAgICAgICAgICAgICAgICAgICAgICAtLSAgIGNvcnJlY3Rpb24gLSBiYXNlIDIKICAgIHwgbiA+IDEgICAgICAgPSAoMi4yNTA2ICwgMC4yODg3ICkgICAgICAgICAgICAgICAgICAgICAgICAgIC0tICAgIChkaXN0LHdpZHRoKQogICAgfCBvdGhlcndpc2UgID0gKDIuMjUwNiAsIDAuNTc3MSApICAKIApudGhIYW0gbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIG46IDEtYmFzZWQgMSwyLDMuLi4KICB8IHcgPj0gMSAgID0gZXJyb3IgJCAiQnJlYWNoIG9mIGNvbnRyYWN0ICh3IDwgMSk6ICAgIiArKyBzaG93ICh3KQogIHwgbSA8ICAwICAgPSBlcnJvciAkICJUaGUgZXN0aW1hdGUgaXMgdG9vIHNtYWxsOiAgICAiICsrIHNob3cgKGMsbikKICB8IG0gPj0gbmIgID0gZXJyb3IgJCAiR2VuZXJhdGVkIGJhbmQgaXMgdG9vIG5hcnJvdzogIiArKyBzaG93IChtLG5iKQogIHwgVHJ1ZSAgICAgPSByZXMgCiB3aGVyZSAKICAoZCx3KSAgID0gcm5ndmFsMiBuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAtLSBjb3JyZWN0aW9uIGRpc3QsIHdpZHRoCiAgaGkgICAgICA9IGVzdHZhbDIgbiAtIGQgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLS0gICBoaSA+IGxvZ3ZhbDIgPiBoaS13IAogIChjLGIpICAgPSBmb2xkbF8gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIHRvdGFsIGNvdW50LCB0aGUgYmFuZAogICAgICAgICAgICAgKFwoIWMsIWIpIChpLHQpIC0+IGNhc2UgdCBvZiB7IFtdIC0+IChpK2MsYikKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgOyBbeF0gLT4gKGkrYyx4OmIpfSkKICAgICAgICAgICAgICgwLFtdKQogICAgICAgICAgICAgWyAoIGkrMSwgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIHRvdGFsIHRyaXBsZXMgdy8gdGhpcyAoaixrKQogICAgICAgICAgICAgICAgIFsgKHIsKGksaixrKSkgfCBmcmFjIDwgdyBdICkgICAgICAgICAgICAgICAgIC0tIHN0b3JlIGl0LCBpZiBpbnNpZGUgYmFuZAogICAgICAgICAgICAgICB8IGsgPC0gWyAwIC4uIGZsb29yICggaGkgICAvbGI1KSBdLCAgbGV0IHAgPSBmcm9tSW50ZWdyYWwgaypsYjUsICAgIAogICAgICAgICAgICAgICAgIGogPC0gWyAwIC4uIGZsb29yICgoaGktcCkvbGIzKSBdLCAgbGV0IHEgPSBmcm9tSW50ZWdyYWwgaipsYjMgKyBwLAogICAgICAgICAgICAgICAgIGxldCAoaSxmcmFjKSA9IHByICAoaGktcSkgOyAgICAgICAgICAgICAgciA9IGhpIC0gZnJhYyAgLS0gciA9IGkgKyBxIAogICAgICAgICAgICAgICBdIHdoZXJlICAgICBwciA9IHByb3BlckZyYWN0aW9uCiAgKG0sbmIpICA9ICggZnJvbUludGVncmFsICQgYyAtIG4sIGxlbmd0aCBiICkgICAgICAgICAgICAgIC0tIG0gMC1iYXNlZCBmcm9tIHRvcCwgfGJhbmR8CiAgKHMscmVzKSA9ICggc29ydEJ5IChmbGlwIGNvbXBhcmUgYG9uYCBmc3QpIGIsIHMhIW0gKSAgICAgIC0tIHNvcnRlZCBkZWNyZWFzaW5nLCByZXN1bHQKICBmb2xkbF8gID0gZm9sZGwn