2

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

![替代文字][1]

哪里q是双，t也是一个双和k一个整数。e代表指数函数。该系数用于步骤中，q并且始终t从k0 开始，直到（该步骤的）所有先前系数的总和达到 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但是当事情变得奇怪时，例如，我开始计算它，我得到的是一系列 0，然后是从某个数字到最后的 NaN。tq = 50.0t = 100.0k = 0 to 100

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

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

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

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

4

2 回答 2

3 
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!]并将其放入数组中，从而避免任何错误。

当然，您可以进一步做多种变化。

于 2010-09-03T19:45:01.503 回答
1

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

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

替代文字

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

于 2010-09-04T00:50:08.763 回答