如何优化此计算

Question

我正在编写一个模型检查器，它依赖于以下算法密集使用的系数的计算：

![替代文字][1]

其中q 是双精度，t 也是双精度，k 是整数。 e 代表指数函数。此系数用于q 和t 不变的步骤，而k 始终从 0 开始，直到（该步骤的）所有先前系数的总和达到 1。

我的第一个实现是字面的：

let rec fact k =
  match k with
    0 | 1 -> 1
    | n -> n * (fact (k - 1))

let coeff q t k = exp(-. q *. t) *. ((q *. t) ** (float k)) /. float (fact k)

当然，这并没有持续多久，因为当k 超过一个小阈值 (15-20) 时计算整个阶乘是不可行的：显然结果开始变得疯狂。所以我通过增量划分重新安排了整个事情：

let rec div_by_fact v d =
  match d with
    1. | 0. -> v
    | d -> div_by_fact (v /. d) (d -. 1.)

let coeff q t k = div_by_fact (exp(-. q *. t) *. ((q *. t) ** (float k))) (float k)

当q 和t 足够“正常”但当事情变得奇怪时，例如q = 50.0 和t = 100.0 我开始从k = 0 to 100 计算它我得到的是一个一系列 0，后跟 NaN，从某个数字开始直到结束。

这当然是由于数字开始太接近于 0 或类似问题的操作造成的。

您对我如何优化公式以便能够在广泛的输入范围内提供足够准确的结果有任何想法吗？

一切都应该是 64 位的（因为我使用的是 OCaml，它默认使用双精度）。也许也有使用 128 位双精度的方法，但我不知道如何。

我正在使用 OCaml，但您可以使用任何您想要的语言提供想法：C、C++、Java 等。我完全使用了所有这些语言。

Answer 1

qt^k/k! = e^[log[qt^k/k!]]
log[qt^k/k!] = log[qt^k] - log[k!] // log[k!] ~ klnk - k  by stirling
             ~ k ln(qt) - (k lnk - k)
             ~ k ln(qt/k) - k

对于较小的 k 值，斯特林近似是不准确的。 但是，由于您似乎在做有限的已知范围，您可以计算 log[k!] 并将其放入数组中，避免任何错误。

当然，您还可以做更多的变化。

Answer 2

这不是一个答案（我相信），而可能只是一个澄清。如果我误解了什么，我会在您发表评论后将其删除。

据我了解，您正在尝试计算 n，例如以下总和等于 1。

您可能会看到它逐渐接近 1，它永远不会等于 1。 如果我误解了你的问题，请纠正我。