fork download
copy 
  1. {-# OPTIONS -O2 -XBangPatterns #-}
  2. import Data.List 
  3. import Data.Function (on)
  4. import Control.Arrow
  5.  
  6. main = let (r,t) = nthHam 10000000000 in print t -- >> print (trival t)
  7.  
  8. lg3 = logBase 2 3;  lg5 = logBase 2 5;  ww = logBase 2 30 / 2
  9. logval (i,j,k)    = fromIntegral i + fromIntegral j*lg3 + fromIntegral k*lg5
  10. trival (i,j,k)    = 2^i * 3^j * 5^k
  11. estval n          = (6*lg3*lg5* fromIntegral n)**(1/3)  -- estimated logval, base 2
  12. rngval n       
  13.     | n > 500000  = (ww - (3/estval n), 6/estval n)     -- the space tweak! (thx, GBG!)
  14.     | n > 500000  = (2.4496 , 0.0076 )                  -- empirical estimation 
  15.     | n > 50000   = (2.4424 , 0.0146 )                  --   correction, base 2
  16.     | n > 500     = (2.3948 , 0.0723 )                  --     (dist,width)
  17.     | n > 1       = (2.2506 , 0.2887 )                  -- around (log $ sqrt 30), 
  18.     | otherwise   = (2.2506 , 0.5771 )                  --   says WP
  19.  
  20. nthHam :: Integer -> (Double, (Int, Int, Int))          -- on 64-bits, use Int!
  21. nthHam n                                                   -- n: 1-based: 1,2,3...
  22.   | w >= 1    = error $ "Breach of contract: (w < 1):  " ++ show w
  23.   | m <  0    = error $ "Not enough triples generated: " ++ show (c,n)
  24.   | m >= nb   = error $ "Generated band is too narrow: " ++ show (m,nb)
  25.   | otherwise = res
  26.  where
  27.   (d,w)   = rngval n                                     -- correction dist, width
  28.   hi      = estval n - d                                 --   hi > logval > hi-w
  29.   (m,nb)  = ( fromIntegral $ c - n, length b )           -- m 0-based from top, |band|
  30.   (s,res) = ( sortBy (flip compare `on` fst) b, s!!m )   -- sorted decreasing, result
  31.   (c,b)   = -- cb hi w
  32.             (sum *** concat) . unzip $
  33.          -- f (0 :: Integer)                             -- total triples count, band
  34.               [ ( fromIntegral i+1,                      -- total triples w/ this (j,k)
  35.                   [ (r,(i,j,k)) | frac < w ] )           -- store it, if inside band
  36.                 | k <- [ 0 .. floor ( hi   /lg5) ],  let p = fromIntegral k*lg5,
  37.                   j <- [ 0 .. floor ((hi-p)/lg3) ],  let q = fromIntegral j*lg3 + p,
  38.                   let (i,frac) = pr  (hi-q) ;            r = hi-frac ] -- r = i + q
  39.          -- f 0 z == (sum $ map fst z, concat $ map snd z)
  40.     where pr = properFraction                            -- pr 12.5 == (12, 0.5)
  41. {-        f !c []          = (c,[])                      -- code as a loop 
  42.           f !c ((c1,b1):r) = let (cr,br) = f (c+c1) r    --   to prevent space leak
  43.                              in case b1 of { [v] -> (cr,v:br) 
  44.                                            ;  _  -> (cr, br) }
  45. -}
  46.  
  47. cb :: Double -> Double -> (Integer, [(Double,(Int,Int,Int))]) -- ~ 20% faster
  48. cb hi w = gk 0 [] 0
  49.   where
  50.   kmax = floor ( hi   /lg5)
  51.   gk c b k | k > kmax  = (c,b)
  52.            | otherwise = gj c b 0
  53.        where
  54.        p    = fromIntegral k*lg5 
  55.        jmax = floor ((hi-p)/lg3)
  56.        gj c b j | j > jmax  = gk c b (k+1)
  57.                 | otherwise = r `seq` c2 `seq` b2 `seq` gj c2 b2 (j+1)
  58.            where
  59.            q = fromIntegral j*lg3 + p
  60.            (i,frac) = properFraction (hi-q)
  61.            r = hi-frac 
  62.            c2 = c + fromIntegral i+1
  63.            b2 = if frac < w then  (r,(i,j,k)):b else b
Success #stdin #stdout 1.93s 7836KB
comments (?)
stdin
copy 
with `cb` (loops):
1T:  10.97s-7.8MB  (1126,16930,40)      t~n^0.67
100B: 2.35s-5.8MB  (10178,1384,279)     t~n^0.73
20B:  0.76s-5.8MB  (4077,927,890)
10B:  0.44s-5.8MB  (2177,8,1659)
1B:   0.02s-4.8MB  (1334,335,404)

with `f` (left fold over a list):
100B: (10178,1384,279) 3.39s-6.8MB vs 2.74s-6.9MB
10B:  (2177,8,1659)    0.64s-5.8MB vs 0.49s-5.9MB   of foldl' version in q3fma
1B:   (1334,335,404)   0.08s-4.8MB vs 0.05s-4.9MB

with HOF (`(sum *** concat) . unzip $` a list):
100B: prog: out of memory (requested 1048576 bytes)  (50B too)
20B: (4077,927,890)  2.44s-7.8MB
10B: (2177,8,1659)   1.91s-7.8MB
1B:  (1334,335,404)  0.34s-7.8MB
stdout
copy 
(2177,8,1659)
https://ideone.com/GXh4P0
language:
Haskell (ghc 8.4.4)
created:
7 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!