R中的模乘和溢出答案

【问题标题】：Modular multiplication and overflow in RR中的模乘和溢出
【发布时间】：2017-08-06 07:01:09
【问题描述】：

我有一个包含大约 1000 个元素的向量，从大约 1 到 10，我希望所有元素的乘积以 10^6 为模。即，我想要类似的东西

x <- as.integer(runif(1000, 1, 10))
prod(x) %% 1000000

但是，由于 prod(x) 的计算结果为 Inf，因此上面的代码只返回 NaN。

在我的实际应用程序中（这不是最小的，所以我不会把它放在这里），我用for 循环完成了这个，比如

x <- as.integer(runif(1000, 1, 10))

for (i in 1:length(x)) {
  if (i == 1) {
    y <- x[i]
  } else {
    y <- y * x[i]
    y <- y %% 1000000
  }
}
y

最初，我只是想寻求帮助，想出一个比for 循环更好的解决方案，假设我应该在 R 中避免这种情况。但是，在写这个问题时，我现在同样感到困惑的是，我似乎无法为 for 循环生成有效的 MWE。上面的代码返回0，但是当我删除as.integer() 时它工作正常。我花了比我愿意承认的时间更长的时间试图弄清楚为什么会这样，但我很难过。

那么两个问题：

为什么上述 MWE 不起作用？（为什么没有as.integer 也能工作？）
在 R 中实现模乘以避免溢出问题的好方法是什么？理想情况下，我会对基本 R 中的解决方案感到好奇，即没有像 gmp 包这样的东西。我试过sapply/lapply，但到目前为止无济于事。

我发现了用其他语言提出的类似问题，但恐怕我很难解读用这些术语表达的答案。

【问题讨论】：

标签： r multiplication modular-arithmetic

【解决方案1】：

让我们看看会发生什么：

set.seed(42)
x <- as.integer(runif(1000, 1, 10))

y <- integer(1000)
for (i in 1:length(x)) {
  if (i == 1) {
    y[i] <- x[i]
  } else {
    y[i] <- y[i-1] * x[i]
    y[i] <- y[i] %% 1000000
  }
}
y[1:30]
# [1]      9     81    243   1944  11664  58320 408240 816480 898880 292160 460800 225600  30400  91200 456000 104000 936000
#[18] 872000 360000 160000 440000 880000 920000 280000 280000 400000 600000 400000      0      0

好的，所以问题出现在第29次迭代：

400000 * x[29]
#[1] 2e+06

为什么浮点数不会发生这种情况？

set.seed(42)
x <- (runif(1000, 1, 10))

y <- integer(1000)
for (i in 1:length(x)) {
  if (i == 1) {
    y[i] <- x[i]
  } else {
    y[i] <- y[i-1] * x[i]
    y[i] <- y[i] %% 1000000
  }
}
y[1:30]
 #[1] 9.233254e+00 8.710356e+01 3.114175e+02 2.638961e+03 1.788083e+04 1.014176e+05 7.737451e+05 7.115236e+05 9.187133e+05
#[10] 7.484847e+05 8.319994e+05 2.167078e+05 3.966498e+04 1.308492e+05 6.752649e+05 3.880944e+05 8.048925e+05 6.559748e+05
#[19] 4.602498e+05 7.812873e+05 1.380611e+05 3.104154e+05 7.312040e+04 6.961072e+05 2.125757e+05 1.963559e+05 8.859247e+05
#[28] 1.076669e+05 5.407815e+05 6.096419e+05

如您所见，您永远不会得到 1e6 的倍数。

请注意，您的循环与矢量化尝试不同，因为您在每次迭代中计算模数。

您可以考虑使用任意精度数字作为替代方案（请参阅包 Rmpfr）。或者使用更智能的算法。

【讨论】：