汉明数和双精度

3
我正在尝试在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

我注意到如果我将数字表示为2、3和5指数的三元组,如data Power = Power { k2 :: !Int, k3 :: !Int, k5 :: !Int },则可以避免(更慢的)任意精度Integer。其中数字被理解为2k2 * 3k3 * 5k5。然后两个Power的比较变成了:
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

因此,粗略地比较 p₁ = 2i₁ * 3j₁ * 5k₁p₂ = 2i₂ * 3j₂ * 5k₂,我们比较 p₁p₂ 的对数,这些对数应该适合于 Double。但实际上,我们做得更好:首先计算它们的 GCD(通过找到相应指数对的 min - 目前仅涉及Int 算术!),将 p₁p₂ 除以 GCD(通过从相应指数中减去min - 也仅涉及Int 算术),然后比较结果的对数。

但是,由于我们使用了 Double,最终会存在精度损失。这就是我提问的原因所在:

  1. 什么时候双精度浮点数的有限精度会影响我?也就是说,如何估计, , 的顺序,使得比较2i * 3j * 5k与具有“类似”指数的数字的结果将变得不可靠?
  2. 我们通过求GCD而进行的除法会降低此任务的指数,这样会如何修改上一个问题的答案?

我做了一个实验,将以这种方式生成的数字与通过任意精度算术生成的数字进行比较,所有Hamming数字都匹配到第1'000'000'000个(这花费了我大约15分钟和600兆字节的RAM来验证)。但这显然不是一个证明。


1
你的问题1是什么?形如2^i•3^j•5^k的最小数x是多少,以便存在另一个数y在该形式下,且x < y,将log x和log y转换为最近的“Double”值得到X和Y,使得Y ≤ X,因此通过比较“Double”中的对数无法区分x和y?问题2类似,只是2、3或5的每个指数在x或y中至多有一个非零?对数使用哪个底数?(底数的影响可能很小,但它可能会产生舍入误差,这可能会影响第一个失败发生的位置。) - Eric Postpischil
或者说,我们在Double中没有直接拥有$x$和$y$的对数,但是我们可以使用Double算术从2、3和5的对数(每个对数乘以指数然后相加)中计算出它们?你是否将2、3和5的对数作为最接近可表示值在Double中(尽管一些数学库可能存在更大的误差,尽管对数比某些超越函数容易计算)? - Eric Postpischil
1
答案是,如果记忆无误的话(但请务必查看罗塞塔代码页面),可能在万亿分之几或更高。您的GCD技巧很好,但不幸的是,总会有一些要比较的三元组没有公共因子,所以最终我猜它并不重要。我记得在这里的某个答案或罗塞塔上提到过这个问题。 - Will Ness
@EricPostpischil 我认为你对q1的表述正是我所想的!关于q2,我认为也是正确的。对数是自然对数,而且我们确实不直接计算x和y的对数,而是依赖于“log(2)”,“log(3)”和“log(5)”。关于最近可表示值,我该如何验证? - 0xd34df00d
2
这个回答直接回答了你的问题。它提到在计算第一万亿个 Hamming 数时使用了14个有效数字。 - Will Ness
显示剩余8条评论
2个回答

2

根据经验,它大约是10万亿个Hamming数或更高。

使用你的GCD技巧在这里不起作用,因为一些相邻的Hamming数之间可能没有共同的因子。

更新:ideone和其他地方尝试在线运行代码,我们得到

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")

0

我猜你可以使用自适应任意精度来计算对数。

如果您选择以2为底数,则log2(2^i)是微不足道的。这消除了1个因素,而且log2比自然对数更容易计算(https://en.wikipedia.org/wiki/Binary_logarithm提供了一个算法,还有Shanks...)。

对于log2(3)和log2(5),您将开发足够的术语来区分两个操作数。我不知道它是否会比在大整数算术中直接求幂3^j和5^k并计算高位导致更多的操作...但是这些可以预先制表达到所需的数字位数。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接