fork download
  1. {-# OPTIONS_GHC -O2 #-}
  2. {-# LANGUAGE BangPatterns #-}
  3. module Main where
  4. import Data.List
  5. import Data.Function
  6.  
  7. main = print $ nthHam 500000000
  8. ln2 = log 2; ln3 = log 3; ln5 = log 5
  9. trival (i,j,k) = 2^i * 3^j * 5^k
  10. logval (i,j,k) = fromIntegral i*ln2 + fromIntegral j*ln3 + fromIntegral k*ln5
  11. estval n = (6*ln2*ln3*ln5* fromIntegral n)**(1/3)
  12. rngval n
  13. | n > 500000 = (1.698 , 0.0050) -- (d,w)
  14. | n > 50000 = (1.693 , 0.0100)
  15. | n > 500 = (1.66 , 0.0500)
  16. | n > 1 = (1.56 , 0.2000)
  17. | otherwise = (1.56 , 0.4000)
  18.  
  19. nthHam n1 -- n1 1-based 2,3...
  20. | w >= ln2 = error $ "Breach of contract: w >= ln2 " ++ show (w,ln2)
  21. | m < 0 = error $ "Not enough triples generated: " ++ show (c,n)
  22. | m >= nh = error $ "Generated band too narrow: " ++ show (m,nh)
  23. | True = ((m,nh) , (\s->and $ zipWith (>) s (tail s)) $ map fst s , res)
  24. where
  25. n = n1 + 1 -- n 1-based 1,2,3...
  26. (s,res) = ( sortBy (flip compare `on` fst) h, s!!m )
  27. (d,w) = rngval n
  28. (w',hi) = ( w/ln2, estval n - d ) -- hi > logval > hi-w
  29. (m,nh) = ( fromInteger $ c - n, length h ) -- m 0-based, from top
  30. {-(c,h) = -- prod (sum,concat) . unzip $ ... -- prod (f,g) (x,y) = (f x,g y)
  31. aux
  32. [ ( i+1, [ (r,(i,j,k)) | frac < w' ])
  33. | k <- [ 0 .. floor ( hi /ln5) ], let p = fromIntegral k*ln5,
  34. j <- [ 0 .. floor ((hi-p)/ln3) ], let q = fromIntegral j*ln3 + p,
  35. let (i,frac) = pr ((hi-q)/ln2) ; r = fromIntegral i*ln2 + q
  36. ] 0 [] where pr = properFraction
  37. aux ((i,[]):t) c b = aux t (c+i) b
  38. aux ((i,[x]):t) c b = aux t (c+i) (x:b)
  39. aux [] c b = (c,b)
  40. -}
  41. (c,h) = auxk [0..floor ( hi /ln5)] 0 []
  42.  
  43. auxk [] c b = (c,b)
  44. auxk (k:ks) c b =
  45. case fromIntegral k*ln5 of { p ->
  46. auxj [0..floor ((hi-p)/ln3)] (k,p,ks) c b }
  47.  
  48. auxj [] (_,_,ks) c b = auxk ks c b
  49. auxj (j:js) st@(k,p,_) c b =
  50. case fromIntegral j*ln3 + p of { q ->
  51. case properFraction ((hi-q)/ln2) of { (i,frac) ->
  52. case {- id $! -} (c+i+1) of { c' ->
  53. {- c' or !c' CHANGES THE WHOLE SHEBANG!
  54. it's a result-producing 0.73s- 4.7MB WITH IT
  55. or stack-overflow 0.94s-74.4MB WITHOUT IT -}
  56. -- try the suggestion by "ehird", in place of the bang pattern,
  57. -- to force the 'c' right here:
  58. -- http://stackoverflow.com/q/9149183/849891 ... Yes, it works!
  59. -- (id $! x == x `seq` x) forces nothing by itself:
  60. -- if (x `seq` y) is forced, then x is forced, and y's the answer
  61. -- and pattern-matching by case {} is what forces stuff
  62. case c' `seq` False of { True -> undefined ; _ ->
  63. if frac < w'
  64. then case (hi-frac*ln2,(i,j,k)):b of { b' ->
  65. {- b' must be built lazily or stack-overflow ensues -}
  66. auxj js st c' b' }
  67. else auxj js st c' b
  68. }}}}
Success #stdin #stdout 0.78s 4740KB
stdin
Standard input is empty
stdout
((2302,4863),True,(1541.7350588298232,(1384,306,153)))