How to optimize this computation

Question

Hello, I'm writing a model-checker which relies on the computation of a coefficient which is used intensively by the algorithms which is the following:

where q is double, t a double too and k an int. e stands for exponential function. This coefficient is used in steps in which q and t don't change while k always starts from 0 until the sum of all previous coefficients (of that step) reaches 1.

My first implementation was a licteral one:

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)

Of course this didn't last so much since computing the whole factorial was just unfeasable when k went over a small threshold (15-20): obviously results started to go crazy. So I rearranged the whole thing by doing incremental divisions:

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)

This version works quite well when q and t are enough 'normal' but when things gets strange, eg q = 50.0 and t = 100.0 and I start to calculate it from k = 0 to 100 what I get is a serie of 0s followed by nans from a certain number until the end.

Of course this is caused by operations with numbers that start to get too near to 0 or similar problems.

Do you have any idea in how I can optimize the formula to be able to give enough accurate results over a wide spread of inputs?

Everything should be already 64 bit (since I'm using OCaml which uses doubles by default). Maybe there is a way to use 128 bit dobules too but I don't know how.

I'm using OCaml but you can provide ideas in whatever language you want: C, C++, Java, etc.. I quite used all of them.

Thanks in advance

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

for small values of k, Stirling approximation is not accurate. however, since you appear to be doing finite known range, you can compute log[k!]and put it in array, avoiding any errors whatsoever.

of course there are multiple variations you can do further.

Answer 2

This is not an answer (I believe), but pehaps just a clarification .If I misunderstood something, I'll delete it after your comment.

As I understand, you are trying to calculate n, such as the following sum is equal to 1.

As you may see it approaches to 1 asymptotically, it'll never be EQUAL to 1. Please correct me if I misunderstood your question.