{-# OPTIONS_GHC -O3 -XStrict #-}
import Data.Word
import Data.List (sortBy)
import Data.Function (on)
nthHam n -- n: 1-based 1,2,3...
| n < 2 = case n of
0 -> error "nthHam: Argument is zero!" _ -> (0, 0, 0) -- trivial case for 1
| m
< 0 = error $ "Not enough triples generated: " ++ show (c
,n
) | m
>= nb
= error $ "Generated band is too narrow: " ++ show (m
,nb
) | otherwise = case res
of (_, tv
) -> tv
-- 2^i * 3^j * 5^k where
lg2est
= (6 * lb3
* lb5
* fromIntegral n
)**(1/3) - lbrt30
-- estimated logval, base 2 (hi,lo) = (lg2est + 1/lg2est, 2 * lg2est - hi) -- hi > log2est > lo
(c
, b
) = let klmt
= floor (hi
/ lb5
) loopk k ck bndk =
if k > klmt then (ck, bndk) else
loopj j cj bndj =
if j > jlmt then loopk (k + 1) cj bndj else
r = hi - frac
if r < lo then bndj
else case (r, (i, j, k)) of
nhd
-> nhd `
seq` nhd : bndj
in ncj `
seq` nbndj `
seq` loopj nj ncj nbndj
in loopj 0 ck bndk
in loopk 0 0 []
-- (s,res) = (b, s!!m)
(s
,res
) = ( sortBy
(flip compare `on`
fst) b
, s
!!m
) -- sorted decreasing, result<
ey0jIE9QVElPTlNfR0hDIC1PMyAtWFN0cmljdCAjLX0KCmltcG9ydCBEYXRhLldvcmQKaW1wb3J0IERhdGEuTGlzdCAoc29ydEJ5KQppbXBvcnQgRGF0YS5GdW5jdGlvbiAob24pCgpudGhIYW0gOjogV29yZDY0IC0+IChJbnQsIEludCwgSW50KQpudGhIYW0gbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLS0gbjogMS1iYXNlZCAxLDIsMy4uLgogIHwgbiA8IDIgICAgID0gY2FzZSBuIG9mCiAgICAgICAgICAgICAgICAgIDAgLT4gZXJyb3IgIm50aEhhbTogIEFyZ3VtZW50IGlzIHplcm8hIgogICAgICAgICAgICAgICAgICBfIC0+ICgwLCAwLCAwKSAgICAgICAgICAgICAgICAgICAgICAgLS0gdHJpdmlhbCBjYXNlIGZvciAxCiAgfCBtIDwgIDAgICAgPSBlcnJvciAkICJOb3QgZW5vdWdoIHRyaXBsZXMgZ2VuZXJhdGVkOiAiICsrIHNob3cgKGMsbikKICB8IG0gPj0gbmIgICA9IGVycm9yICQgIkdlbmVyYXRlZCBiYW5kIGlzIHRvbyBuYXJyb3c6ICIgKysgc2hvdyAobSxuYikKICB8IG90aGVyd2lzZSA9IGNhc2UgcmVzIG9mIChfLCB0dikgLT4gdHYgLS0gMl5pICogM15qICogNV5rCiB3aGVyZQogIGxiMyAgICAgPSBsb2dCYXNlIDIgMzsgbGI1ID0gbG9nQmFzZSAyIDUuMAogIGxicnQzMCAgPSBsb2dCYXNlIDIgJCBzcXJ0IDMwIDo6IERvdWJsZSAtLSBlc3RpbWF0ZSBhZGp1c3RtZW50IGFzIHBlciBXUAogIGxnMmVzdCAgPSAoNiAqIGxiMyAqIGxiNSAqIGZyb21JbnRlZ3JhbCBuKSoqKDEvMykgLSBsYnJ0MzAgLS0gZXN0aW1hdGVkIGxvZ3ZhbCwgYmFzZSAyCiAgKGhpLGxvKSA9IChsZzJlc3QgKyAxL2xnMmVzdCwgMiAqIGxnMmVzdCAtIGhpKSAtLSBoaSA+IGxvZzJlc3QgPiBsbwogIChjLCBiKSAgPSBsZXQga2xtdCA9IGZsb29yIChoaSAvIGxiNSkKICAgICAgICAgICAgICAgIGxvb3BrIGsgY2sgYm5kayA9CiAgICAgICAgICAgICAgICAgIGlmIGsgPiBrbG10IHRoZW4gKGNrLCBibmRrKSBlbHNlCiAgICAgICAgICAgICAgICAgIGxldCBwID0gaGkgLSBmcm9tSW50ZWdyYWwgayAqIGxiNTsgamxtdCA9IGZsb29yIChwIC8gbGIzKQogICAgICAgICAgICAgICAgICAgICAgbG9vcGogaiBjaiBibmRqID0KICAgICAgICAgICAgICAgICAgICAgICAgaWYgaiA+IGpsbXQgdGhlbiBsb29wayAoayArIDEpIGNqIGJuZGogZWxzZQogICAgICAgICAgICAgICAgICAgICAgICBsZXQgcSA9IHAgLSBmcm9tSW50ZWdyYWwgaiAqIGxiMwogICAgICAgICAgICAgICAgICAgICAgICAgICAgKGksIGZyYWMpID0gcHJvcGVyRnJhY3Rpb24gcQogICAgICAgICAgICAgICAgICAgICAgICAgICAgbmogPSBqICsgMTsgbmNqID0gY2ogKyBmcm9tSW50ZWdyYWwgaSArIDEKICAgICAgICAgICAgICAgICAgICAgICAgICAgIHIgPSBoaSAtIGZyYWMKICAgICAgICAgICAgICAgICAgICAgICAgICAgIG5ibmRqID0gaSBgc2VxYCBibmRqIGBzZXFgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGlmIHIgPCBsbyB0aGVuIGJuZGoKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZWxzZSBjYXNlIChyLCAoaSwgaiwgaykpIG9mCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuaGQgLT4gbmhkIGBzZXFgIG5oZCA6IGJuZGoKICAgICAgICAgICAgICAgICAgICAgICAgaW4gbmNqIGBzZXFgIG5ibmRqIGBzZXFgIGxvb3BqIG5qIG5jaiBuYm5kagogICAgICAgICAgICAgICAgICBpbiBsb29waiAwIGNrIGJuZGsKICAgICAgICAgICAgaW4gbG9vcGsgMCAwIFtdCiAgKG0sbmIpICA9ICggZnJvbUludGVncmFsICQgYyAtIG4sIGxlbmd0aCBiICkgICAgICAgICAtLSBtIDAtYmFzZWQgZnJvbSB0b3AsIHxiYW5kfAotLSAgKHMscmVzKSA9IChiLCBzISFtKQogIChzLHJlcykgPSAoIHNvcnRCeSAoZmxpcCBjb21wYXJlIGBvbmAgZnN0KSBiLCBzISFtICkgLS0gc29ydGVkIGRlY3JlYXNpbmcsIHJlc3VsdDwKCm1haW4gPSBwdXRTdHJMbiAkIHNob3cgJCBudGhIYW0gMTAwMDAwMDAwMDAwMA==