{-# OPTIONS -O2 #-}
-- directly find n-th Hamming number, 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
import Data.List -- cf. ideone.com/q3fma
import Data.Function
main
= let (r
,t
) = nthHam
1000000000000 in print t
-- >> print (trival t)
trival (i,j,k) = 2^i * 3^j * 5^k
estval n
| n > 500000 = (v - lb30_2 + (3/v), 6/v) -- the space tweak! (thx, GBG!)
| n > 500000 = (v - 2.4496 , 0.0076 ) -- empirical estimation
| n > 50000 = (v - 2.4424 , 0.0146 ) -- correction, base 2
| n > 500 = (v - 2.3948 , 0.0723 ) -- (dist,width)
| n > 1 = (v - 2.2506 , 0.2887 ) -- around (log $ sqrt 30),
| otherwise = (v
- 2.2506 , 0.5771 ) -- says WP where v
= (6*lb3
*lb5
* fromIntegral n
)**(1/3) -- estimated logval, base 2
nthHam n -- n: 1-based: 1,2,3...
| n
<= 0 = error $ "n is 1--based: must be n > 0: " ++ show n
| 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
) where
(hi,w) = estval n -- hi > logval > hi-w
nb
= length b
-- length of the band {- (Sum c,b) = fold [ (Sum (fromIntegral i+1), ... ] -}
(c
,b
) = foldl
_ (\
(c
,b
) (i
,t
)-> let c2
=c
+i
in c2`
seq`
-- ( total count, the band ) case t of []-> (c2,b);[v]->(c2,v:b) ) (0,[])
-- prod (sum,concat) . unzip $
[ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band
let (i,frac) = pr (hi-q) ; r = hi - frac -- r = i + q
prod (f,g) (x,y) = (f x,g y) -- (f *** g)
ey0jIE9QVElPTlMgLU8yICMtfQoKLS0gZGlyZWN0bHkgZmluZCBuLXRoIEhhbW1pbmcgbnVtYmVyLCBpbiB+IE8obl57Mi8zfSkgdGltZQotLSBiYXNlZCBvbiAidG9wIGJhbmQiIGlkZWEgYnkgTG91aXMgS2xhdWRlciwgZnJvbSBEREogZGlzY3Vzc2lvbgotLSBieSBXaWxsIE5lc3MsIG9yaWdpbmFsIHBvc3Q6IGRyZG9iYnMuY29tL2Jsb2dzL2FyY2hpdGVjdHVyZS1hbmQtZGVzaWduLzIyODcwMDUzOAoKaW1wb3J0IERhdGEuTGlzdCAgICAgICAgICAgLS0gY2YuIGlkZW9uZS5jb20vcTNmbWEKaW1wb3J0IERhdGEuRnVuY3Rpb24gIAoKbWFpbiA9IGxldCAocix0KSA9IG50aEhhbSAxMDAwMDAwMDAwMDAwIGluIHByaW50IHQgLS0gPj4gcHJpbnQgKHRyaXZhbCB0KQoKbGIzID0gbG9nQmFzZSAyIDM7ICBsYjUgPSBsb2dCYXNlIDIgNTsgIGxiMzBfMiA9IGxvZ0Jhc2UgMiAzMCAvIDIKdHJpdmFsIChpLGosaykgICAgPSAyXmkgKiAzXmogKiA1XmsKZXN0dmFsIG4gICAgICAgCiAgICB8IG4gPiA1MDAwMDAgID0gKHYgLSBsYjMwXzIgKyAoMy92KSwgNi92KSAgICAgICAgICAgICAtLSB0aGUgc3BhY2UgdHdlYWshICh0aHgsIEdCRyEpCiAgICB8IG4gPiA1MDAwMDAgID0gKHYgLSAyLjQ0OTYgLCAwLjAwNzYgKSAgICAgICAgICAgICAgICAtLSBlbXBpcmljYWwgZXN0aW1hdGlvbiAKICAgIHwgbiA+IDUwMDAwICAgPSAodiAtIDIuNDQyNCAsIDAuMDE0NiApICAgICAgICAgICAgICAgIC0tICAgY29ycmVjdGlvbiwgYmFzZSAyCiAgICB8IG4gPiA1MDAgICAgID0gKHYgLSAyLjM5NDggLCAwLjA3MjMgKSAgICAgICAgICAgICAgICAtLSAgICAgKGRpc3Qsd2lkdGgpCiAgICB8IG4gPiAxICAgICAgID0gKHYgLSAyLjI1MDYgLCAwLjI4ODcgKSAgICAgICAgICAgICAgICAtLSBhcm91bmQgKGxvZyAkIHNxcnQgMzApLCAKICAgIHwgb3RoZXJ3aXNlICAgPSAodiAtIDIuMjUwNiAsIDAuNTc3MSApICAgICAgICAgICAgICAgIC0tICAgc2F5cyBXUAogIHdoZXJlIHYgPSAoNipsYjMqbGI1KiBmcm9tSW50ZWdyYWwgbikqKigxLzMpICAgICAgICAgICAgLS0gZXN0aW1hdGVkIGxvZ3ZhbCwgYmFzZSAyCgpudGhIYW0gOjogSW50IC0+IChEb3VibGUsIChJbnQsIEludCwgSW50KSkgICAgICAgICAgICAgICAgLS0gKCA2NGJpdDogdXNlIEludCEhISAgICAgTkIhICkKbnRoSGFtIG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIG46IDEtYmFzZWQ6IDEsMiwzLi4uCiAgfCBuIDw9IDAgICAgPSBlcnJvciAkICJuIGlzIDEtLWJhc2VkOiBtdXN0IGJlIG4gPiAwOiAiICsrIHNob3cgbgogIHwgdyA+PSAxICAgID0gZXJyb3IgJCAiQnJlYWNoIG9mIGNvbnRyYWN0OiAodyA8IDEpOiAgIiArKyBzaG93IHcKICB8IG0gPCAgMCAgICA9IGVycm9yICQgIk5vdCBlbm91Z2ggdHJpcGxlcyBnZW5lcmF0ZWQ6ICIgKysgc2hvdyAoYyxuKQogIHwgbSA+PSBuYiAgID0gZXJyb3IgJCAiR2VuZXJhdGVkIGJhbmQgaXMgdG9vIG5hcnJvdzogIiArKyBzaG93IChtLG5iKQogIHwgb3RoZXJ3aXNlID0gc29ydEJ5IChmbGlwIGNvbXBhcmUgYG9uYCBmc3QpIGIgISEgbSAgICAgLS0gbS10aCBmcm9tIHRvcCBpbiBzb3J0ZWQgYmFuZAogd2hlcmUgIAogIChoaSx3KSA9IGVzdHZhbCBuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLS0gICBoaSA+IGxvZ3ZhbCA+IGhpLXcKICBtICAgICAgPSBmcm9tSW50ZWdyYWwgKGMgLSBuKSAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIHRhcmdldCBpbmRleCwgZnJvbSB0b3AKICBuYiAgICAgPSBsZW5ndGggYiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIGxlbmd0aCBvZiB0aGUgYmFuZAogIHstIChTdW0gYyxiKSA9IGZvbGQgWyAoU3VtIChmcm9tSW50ZWdyYWwgaSsxKSwgLi4uIF0gLX0gICAgICAKICAoYyxiKSAgPSBmb2xkbF8gKFwoYyxiKSAoaSx0KS0+IGxldCBjMj1jK2kgaW4gYzJgc2VxYCAgIC0tICggdG90YWwgY291bnQsIHRoZSBiYW5kICkKICAgICAgICAgICAgICAgICAgICAgY2FzZSB0IG9mIFtdLT4gKGMyLGIpO1t2XS0+KGMyLHY6YikgKSAoMCxbXSkKICAgICAgICAgICAtLSBwcm9kIChzdW0sY29uY2F0KSAuIHVuemlwICQgCiAgICAgICAgICAgIFsgKCBmcm9tSW50ZWdyYWwgaSsxLCAgICAgICAgICAgICAgICAgICAgICAgICAtLSB0b3RhbCB0cmlwbGVzIHcvIHRoaXMgKGosaykKICAgICAgICAgICAgICAgIFsgKHIsKGksaixrKSkgfCBmcmFjIDwgdyBdICkgICAgICAgICAgICAgIC0tIHN0b3JlIGl0LCBpZiBpbnNpZGUgYmFuZAogICAgICAgICAgICAgIHwgayA8LSBbIDAgLi4gZmxvb3IgKCBoaSAgIC9sYjUpIF0sICBsZXQgcCA9IGZyb21JbnRlZ3JhbCBrKmxiNSwgICAgCiAgICAgICAgICAgICAgICBqIDwtIFsgMCAuLiBmbG9vciAoKGhpLXApL2xiMykgXSwgIGxldCBxID0gZnJvbUludGVncmFsIGoqbGIzICsgcCwKICAgICAgICAgICAgICAgIGxldCAoaSxmcmFjKSA9IHByICAoaGktcSkgOyAgICAgICAgICAgIHIgPSBoaSAtIGZyYWMgIC0tIHIgPSBpICsgcSAKICAgICAgICAgICAgXSB3aGVyZSAgcHIgICAgICA9IHByb3BlckZyYWN0aW9uICAgICAgICAgICAgIC0tIHByIDEuMjQgPT4gKDEsMC4yNCkKICAgICAgICAgICAgICAgICAgICAgcHJvZCAoZixnKSAoeCx5KSA9IChmIHgsZyB5KSAgICAgICAgIC0tIChmICoqKiBnKQogICAgICAgICAgICAgICAgICAgICBmb2xkbF8gID0gZm9sZGwn
IAoyMDIwLTAzLTIzOiAKbm93IHVzaW5nIEludCBpbiBudGhIYW0gOjogSW50IC0mZ3Q7IChEb3VibGUsIChJbnQsIEludCwgSW50KSkgCmluc3RlYWQgb2YgdGhlIHByZXZpb3VzIEludGVnZXI6CiAgICBudGhIYW0gOjogSW50ZWdlciAtJmd0OyAoRG91YmxlLCAoSW50LCBJbnQsIEludCkpIAoKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBPY3QtMjAxNzogCjFUOiAgICAgICAgICAgICAgKDExMjYsMTY5MzAsNDApICAgICAxLjc5cyAgICBuXjAuNjcKMTAwQjogICAgICAgICAgICAoMTAxNzgsMTM4NCwyNzkpICAgIDAuMzhzCjIwQjogMC44MHMtNS44TUIgKDQwNzcsOTI3LDg5MCkgICAgICAwLjE0cwoxMEI6IDAuNDVzLTUuOE1CICgyMTc3LDgsMTY1OSkgICAgICAgMC4wOHMKMUI6ICAwLjA1cy00LjhNQiAoMTMzNCwzMzUsNDA0KQ==
2020-03-23:
now using Int in nthHam :: Int -> (Double, (Int, Int, Int))
instead of the previous Integer:
nthHam :: Integer -> (Double, (Int, Int, Int))
Oct-2017:
1T: (1126,16930,40) 1.79s n^0.67
100B: (10178,1384,279) 0.38s
20B: 0.80s-5.8MB (4077,927,890) 0.14s
10B: 0.45s-5.8MB (2177,8,1659) 0.08s
1B: 0.05s-4.8MB (1334,335,404)