{-# OPTIONS_GHC -O2 #-}
{-# LANGUAGE BangPatterns #-}
module Main where
import Data.List
import Data.Function

main = print $ nthHam 500000000
ln2 = log 2;  ln3 = log 3;  ln5 = log 5  
trival (i,j,k)    =  2^i  *  3^j  *  5^k          
logval (i,j,k)    = fromIntegral i*ln2 + fromIntegral j*ln3 + fromIntegral k*ln5
estval n          = (6*ln2*ln3*ln5* fromIntegral n)**(1/3)
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)  
    | otherwise  = (1.56  , 0.4000)  

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)
  | True     = ((m,nh) , (\s->and $ zipWith (>) s (tail s)) $ map fst s , res)
 where 
  n       = n1 + 1                                      -- n 1-based 1,2,3...
  (s,res) = ( sortBy (flip compare `on` fst) h, s!!m )
  (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' ])    
               | k <- [ 0 .. floor ( hi   /ln5) ],  let p = fromIntegral k*ln5,    
                 j <- [ 0 .. floor ((hi-p)/ln3) ],  let q = fromIntegral j*ln3 + p,
                 let (i,frac) = pr ((hi-q)/ln2) ;       r = fromIntegral i*ln2 + q 
               ] 0 []    where  pr = properFraction 
  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 = 
    case fromIntegral k*ln5 of { p ->
          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 fromIntegral j*ln3 + p of { q ->
     case properFraction ((hi-q)/ln2) of { (i,frac) ->
      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
          }}}}