{-# OPTIONS_GHC -O2 #-}
{-# LANGUAGE BangPatterns #-}
module Main where
import Data.List
import Data.Function
main
= print $ nthHam
500000000trival (i,j,k) = 2^i * 3^j * 5^k
rngval n
| n > 500000 = (1.698 , 0.0050) -- (d,w)
| n > 50000 = (1.693 , 0.0100)
| n > 500 = (1.66 , 0.0500)
| n > 1 = (1.56 , 0.2000)
nthHam n1 -- n1 1-based 2,3...
| w
>= ln2
= error $ "Breach of contract: w >= ln2 " ++ show (w
,ln2
) | m
< 0 = error $ "Not enough triples generated: " ++ show (c
,n
) | m
>= nh
= error $ "Generated band too narrow: " ++ show (m
,nh
) where
n = n1 + 1 -- n 1-based 1,2,3...
(d,w) = rngval n
(w',hi) = ( w/ln2, estval n - d ) -- hi > logval > hi-w
(m,nh) = ( fromInteger $ c - n, length h ) -- m 0-based, from top
{-(c,h) = -- prod (sum,concat) . unzip $ ... -- prod (f,g) (x,y) = (f x,g y)
aux
[ ( i+1, [ (r,(i,j,k)) | frac < w' ])
aux ((i,[]):t) c b = aux t (c+i) b
aux ((i,[x]):t) c b = aux t (c+i) (x:b)
aux [] c b = (c,b)
-}
(c
,h
) = auxk
[0..floor ( hi
/ln5
)] 0 []
auxk [] c b = (c,b)
auxk (k:ks) c b =
auxj
[0..floor ((hi
-p
)/ln3
)] (k
,p
,ks
) c b
}
auxj [] (_,_,ks) c b = auxk ks c b
auxj (j:js) st@(k,p,_) c b =
case {- id $! -} (c+i+1) of { c' ->
{- c' or !c
' CHANGES THE WHOLE SHEBANG! it's a result-producing 0.73s- 4.7MB WITH IT
or stack
-overflow 0
.94s
-74
.4MB WITHOUT IT
-} -- try the suggestion by "ehird", in place of the bang pattern,
-- to force the 'c' right here:
-- http://stackoverflow.com/q/9149183/849891 ... Yes, it works!
-- (id $! x == x `seq` x) forces nothing by itself:
-- if (x `seq` y) is forced, then x is forced, and y's the answer
-- and pattern-matching by case {} is what forces stuff
case c' `seq` False of { True -> undefined ; _ ->
if frac < w'
then case (hi-frac*ln2,(i,j,k)):b of { b' ->
{- b' must be built lazily
or stack
-overflow ensues
-} auxj js st c' b' }
else auxj js st c' b
}}}}
ey0jIE9QVElPTlNfR0hDIC1PMiAjLX0Key0jIExBTkdVQUdFIEJhbmdQYXR0ZXJucyAjLX0KbW9kdWxlIE1haW4gd2hlcmUKaW1wb3J0IERhdGEuTGlzdAppbXBvcnQgRGF0YS5GdW5jdGlvbgoKbWFpbiA9IHByaW50ICQgbnRoSGFtIDUwMDAwMDAwMApsbjIgPSBsb2cgMjsgIGxuMyA9IGxvZyAzOyAgbG41ID0gbG9nIDUgIAp0cml2YWwgKGksaixrKSAgICA9ICAyXmkgICogIDNeaiAgKiAgNV5rICAgICAgICAgIApsb2d2YWwgKGksaixrKSAgICA9IGZyb21JbnRlZ3JhbCBpKmxuMiArIGZyb21JbnRlZ3JhbCBqKmxuMyArIGZyb21JbnRlZ3JhbCBrKmxuNQplc3R2YWwgbiAgICAgICAgICA9ICg2KmxuMipsbjMqbG41KiBmcm9tSW50ZWdyYWwgbikqKigxLzMpCnJuZ3ZhbCBuICAgICAgIAogICAgfCBuID4gNTAwMDAwICA9ICgxLjY5OCAsIDAuMDA1MCkgICAgICAgICAgICAgICAgICAgLS0gKGQsdykKICAgIHwgbiA+IDUwMDAwICAgPSAoMS42OTMgLCAwLjAxMDApICAKICAgIHwgbiA+IDUwMCAgICAgPSAoMS42NiAgLCAwLjA1MDApICAKICAgIHwgbiA+IDEgICAgICAgPSAoMS41NiAgLCAwLjIwMDApICAKICAgIHwgb3RoZXJ3aXNlICA9ICgxLjU2ICAsIDAuNDAwMCkgIAoKbnRoSGFtIG4xICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAtLSBuMSAxLWJhc2VkIDIsMy4uLgogIHwgdyA+PSBsbjIgPSBlcnJvciAkICJCcmVhY2ggb2YgY29udHJhY3Q6IHcgPj0gbG4yICIgICsrIHNob3cgKHcsbG4yKQogIHwgbSA8ICAwICAgPSBlcnJvciAkICJOb3QgZW5vdWdoIHRyaXBsZXMgZ2VuZXJhdGVkOiAiICsrIHNob3cgKGMsbikKICB8IG0gPj0gbmggID0gZXJyb3IgJCAiR2VuZXJhdGVkIGJhbmQgdG9vIG5hcnJvdzogIiAgICArKyBzaG93IChtLG5oKQogIHwgVHJ1ZSAgICAgPSAoKG0sbmgpICwgKFxzLT5hbmQgJCB6aXBXaXRoICg+KSBzICh0YWlsIHMpKSAkIG1hcCBmc3QgcyAsIHJlcykKIHdoZXJlIAogIG4gICAgICAgPSBuMSArIDEgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0tIG4gMS1iYXNlZCAxLDIsMy4uLgogIChzLHJlcykgPSAoIHNvcnRCeSAoZmxpcCBjb21wYXJlIGBvbmAgZnN0KSBoLCBzISFtICkKICAoZCx3KSAgID0gcm5ndmFsIG4KICAodycsaGkpID0gKCB3L2xuMiwgZXN0dmFsIG4gLSBkICkgICAgICAgICAgICAgICAgICAgICAtLSBoaSA+IGxvZ3ZhbCA+IGhpLXcKICAobSxuaCkgID0gKCBmcm9tSW50ZWdlciAkIGMgLSBuLCBsZW5ndGggaCApICAgICAgICAgLS0gbSAwLWJhc2VkLCBmcm9tIHRvcCAKICB7LShjLGgpICAgPSAtLSBwcm9kIChzdW0sY29uY2F0KSAuIHVuemlwICQgLi4uICAgIC0tIHByb2QgKGYsZykgKHgseSkgPSAoZiB4LGcgeSkKICAgICAgICAgICAgICBhdXgKICAgICAgICAgICAgIFsgKCBpKzEsIFsgKHIsKGksaixrKSkgfCBmcmFjIDwgdycgXSkgICAgCiAgICAgICAgICAgICAgIHwgayA8LSBbIDAgLi4gZmxvb3IgKCBoaSAgIC9sbjUpIF0sICBsZXQgcCA9IGZyb21JbnRlZ3JhbCBrKmxuNSwgICAgCiAgICAgICAgICAgICAgICAgaiA8LSBbIDAgLi4gZmxvb3IgKChoaS1wKS9sbjMpIF0sICBsZXQgcSA9IGZyb21JbnRlZ3JhbCBqKmxuMyArIHAsCiAgICAgICAgICAgICAgICAgbGV0IChpLGZyYWMpID0gcHIgKChoaS1xKS9sbjIpIDsgICAgICAgciA9IGZyb21JbnRlZ3JhbCBpKmxuMiArIHEgCiAgICAgICAgICAgICAgIF0gMCBbXSAgICB3aGVyZSAgcHIgPSBwcm9wZXJGcmFjdGlvbiAKICBhdXggKChpLFtdKTp0KSAgYyBiID0gYXV4IHQgKGMraSkgYiAKICBhdXggKChpLFt4XSk6dCkgYyBiID0gYXV4IHQgKGMraSkgKHg6YikKICBhdXggW10gICAgICAgICAgYyBiID0gKGMsYikKICAtfQogIChjLGgpID0gYXV4ayBbMC4uZmxvb3IgKCBoaSAgL2xuNSldIDAgW10KCiAgYXV4ayBbXSBjIGIgICAgID0gKGMsYikKICBhdXhrIChrOmtzKSBjIGIgPSAKICAgIGNhc2UgZnJvbUludGVncmFsIGsqbG41IG9mIHsgcCAtPgogICAgICAgICAgYXV4aiBbMC4uZmxvb3IgKChoaS1wKS9sbjMpXSAoayxwLGtzKSBjIGIgfQogIAogIGF1eGogW10gICAgIChfLF8sa3MpICAgYyBiID0gYXV4ayBrcyBjIGIKICBhdXhqIChqOmpzKSBzdEAoayxwLF8pIGMgYiA9CiAgICBjYXNlIGZyb21JbnRlZ3JhbCBqKmxuMyArIHAgb2YgeyBxIC0+CiAgICAgY2FzZSBwcm9wZXJGcmFjdGlvbiAoKGhpLXEpL2xuMikgb2YgeyAoaSxmcmFjKSAtPgogICAgICBjYXNlIHstIGlkICQhIC19IChjK2krMSkgb2YgeyBjJyAtPiAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgey0gYycgb3IgIWMnIENIQU5HRVMgVEhFIFdIT0xFIFNIRUJBTkchCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGl0J3MgIGEgcmVzdWx0LXByb2R1Y2luZyAgMC43M3MtIDQuN01CICAgV0lUSCAgSVQKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgb3IgICAgc3RhY2stb3ZlcmZsb3cgICAgICAwLjk0cy03NC40TUIgICBXSVRIT1VUIElUIC19CiAgICAgICAtLSB0cnkgdGhlIHN1Z2dlc3Rpb24gYnkgImVoaXJkIiwgaW4gcGxhY2Ugb2YgdGhlIGJhbmcgcGF0dGVybiwgCiAgICAgICAtLSAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRvIGZvcmNlIHRoZSAnYycgcmlnaHQgaGVyZToKICAgICAgIC0tICAgaHR0cDovL3N0YWNrb3ZlcmZsb3cuY29tL3EvOTE0OTE4My84NDk4OTEgLi4uIFllcywgaXQgd29ya3MhCiAgICAgICAtLSAgICAgICAgICAgKGlkICQhIHggPT0geCBgc2VxYCB4KSBmb3JjZXMgbm90aGluZyBieSBpdHNlbGY6CiAgICAgICAtLSAgICAgICAgICAgaWYgKHggYHNlcWAgeSkgaXMgZm9yY2VkLCB0aGVuIHggaXMgZm9yY2VkLCBhbmQgeSdzIHRoZSBhbnN3ZXIKICAgICAgIC0tICAgICAgICAgICAgIGFuZCBwYXR0ZXJuLW1hdGNoaW5nIGJ5IGNhc2Uge30gaXMgd2hhdCBmb3JjZXMgc3R1ZmYKICAgICAgIGNhc2UgYycgYHNlcWAgRmFsc2Ugb2YgeyBUcnVlIC0+IHVuZGVmaW5lZCA7IF8gLT4gIAogICAgICAgIGlmIGZyYWMgPCB3JwogICAgICAgICAgdGhlbiBjYXNlIChoaS1mcmFjKmxuMiwoaSxqLGspKTpiIG9mIHsgYicgLT4KICAgICAgICAgICAgICAgICAgICB7LSBiJyBtdXN0IGJlIGJ1aWx0IGxhemlseSBvciBzdGFjay1vdmVyZmxvdyBlbnN1ZXMgLX0KICAgICAgICAgICAgICAgIGF1eGoganMgc3QgYycgYicgfQogICAgICAgICAgZWxzZSAgYXV4aiBqcyBzdCBjJyBiCiAgICAgICAgICB9fX19