汉明数和双精度

我正在玩生成汉明数 http://en.wikipedia.org/wiki/Regular_number在 Haskell 中，尝试改进明显的功能（请原谅函数的命名）

mergeUniq :: Ord a => [a] -> [a] -> [a]
mergeUniq (x:xs) (y:ys) = case x `compare` y of
                               EQ -> x : mergeUniq xs ys
                               LT -> x : mergeUniq xs (y:ys)
                               GT -> y : mergeUniq (x:xs) ys

powers :: [Integer]
powers = 1 : expand 2 `mergeUniq` expand 3 `mergeUniq` expand 5
  where
    expand factor = (factor *) <$> powers

I noticed that I can avoid the (slower) arbitrary precision Integer if I represent the numbers as the triple of the 2-, 3- and 5-exponents like data Power = Power { k2 :: !Int, k3 :: !Int, k5 :: !Int }, where the number is understood to be 2k2 * 3k3 * 5k5. The comparison of two Powers then becomes

instance Ord Power where
  p1 `compare` p2 = toComp (p1 `divP` gcdP) `compare` toComp (p2 `divP` gcdP)
    where
    divP p1 p2 = Power { k2 = k2 p1 - k2 p2, k3 = k3 p1 - k3 p2, k5 = k5 p1 - k5 p2 }
    gcdP = Power { k2 = min (k2 p1) (k2 p2), k3 = min (k3 p1) (k3 p2), k5 = min (k5 p1) (k5 p2) }
    toComp Power { .. } = fromIntegral k2 * log 2 + fromIntegral k3 * log 3 + fromIntegral k5 * log 5

So, very roughly speaking, to compare p₁ = 2i₁ * 3j₁ * 5k₁ and p₂ = 2i₂ * 3j₂ * 5k₂ we compare the logarithms of p₁ and p₂, which presumably fit Double. But actually we do even better: we first compute their GCD (by finding the mins of the corresponding exponents pairs — only Int arithmetic so far!), divide p₁ and p₂ by the GCD (by subtracting the mins from the corresponding exponents — also only Int arithmetic), and compare the logarithms of the results.

但是，考虑到我们经历了Doubles，最终会损失精度。这就是我提出问题的基础：

1. When will the finite precision of Doubles bite me? That is, how to estimate the order of i, j, k for which the results of comparisons of 2i * 3j * 5k with numbers with "similar" exponents will become unreliable?
2. 我们除以 GCD（这可能会大大降低此任务的指数）这一事实如何修改前一个问题的答案？

我做了一个实验，将这种方式生成的数字与通过任意精度算术生成的数字进行比较，直到第 1'000'000'000 个的所有汉明数都完全匹配（这花了我大约 15 分钟和 600 兆 RAM核实）。但这显然不是证据。

经验上 https://stackoverflow.com/a/12041774/849891，它高于大约 10 万亿分之一汉明数，或更高。

使用你漂亮的 GCD 技巧在这里对我们没有帮助，因为一些相邻的汉明数之间肯定没有公因数。

update:在线尝试关于ideone https://ideone.com/gHPmTo在其他地方，我们得到

4T  5.81s 22.2MB  -- 16 digits used.... still good
                  --  (as evidenced by the `True` below), but really pushing it.
((True,44531.6794,7.275957614183426e-11),(16348,16503,873),"2.3509E+13405")
-- isTruly  max        min logval           nth-Hamming       approx.
--  Sorted   logval      difference          as i,j,k          value
--            in band      in band                             in decimal
10T   11.13s 26.4MB
((True,60439.6639,7.275957614183426e-11),(18187,23771,1971),"1.4182E+18194")
13T   14.44s 30.4MB    ...still good
((True,65963.6432,5.820766091346741e-11),(28648,21308,1526),"1.0845E+19857")

---- same code on tio:
10T   16.77s
35T   38.84s 
((True,91766.4800,5.820766091346741e-11),(13824,2133,32112),"2.9045E+27624")
70T   59.57s
((True,115619.1575,5.820766091346741e-11),(13125,13687,34799),"6.8310E+34804")

---- on home machine:
100T: 368.13s
((True,130216.1408,5.820766091346741e-11),(88324,876,17444),"9.2111E+39198")

140T: 466.69s
((True,145671.6480,5.820766091346741e-11),(9918,24002,42082),"3.4322E+43851")

170T: 383.26s         ---FAULTY---
((False,155411.2501,0.0),(77201,27980,14584),"2.80508E+46783")