Comment optimiser ce calcul

Question

J'écris un modèle vérificateur qui repose sur le calcul d'un coefficient qui est utilisé intensivement par les algorithmes qui est le suivant:

[alt texte] [1]

où q est double, t un double aussi et k un int. e est synonyme de fonction exponentielle. Ce coefficient est utilisé dans les étapes où q et t ne changent pas pendant k commence toujours de 0 jusqu'à ce que la somme de tous les coefficients précédents (de cette étape) atteint 1.

Ma première mise en œuvre était un littéral:

Bien sûr, cela n'a pas duré autant puisque le calcul de la factorielle entière était tout simplement irréalisable quand k dépassait un petit seuil (15-20): évidemment, les résultats ont commencé à devenir fous. Donc, je réarrangé le tout en faisant des divisions supplémentaires:

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) 

Cette version fonctionne très bien quand q et t sont assez « normal », mais quand les choses devient étrange, par exemple q = 50.0 et t = 100.0 et je commence à calculer à partir k = 0 to 100 ce que je reçois est une série de 0 suivi de NaNs d'un certain nombre jusqu'à la fin.

Bien sûr, cela est dû à des opérations dont les nombres commencent à être trop proches de zéro ou à des problèmes similaires. Avez-vous une idée de la façon dont je peux optimiser la formule pour pouvoir donner des résultats suffisamment précis sur une large gamme d'entrées?

Tout devrait déjà être en 64 bits (puisque j'utilise OCaml qui utilise les doubles par défaut). Peut-être existe-t-il un moyen d'utiliser des 128 bits doubles, mais je ne sais pas comment. J'utilise OCaml mais vous pouvez fournir des idées dans la langue que vous voulez: C, C++, Java, etc. Je les ai toutes utilisées.

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 

Pour de petites valeurs de k, l'approximation de Stirling n'est pas précise. Cependant, comme vous semblez avoir une portée connue limitée, vous pouvez calculer log[k!] et le mettre en tableau, en évitant toute erreur.

Bien sûr, il existe plusieurs variantes que vous pouvez faire plus loin.

Answer 2

Ce n'est pas une réponse (je crois), mais peut-être juste une clarification. Si j'ai mal compris quelque chose, je vais le supprimer après votre commentaire.

Si je comprends bien, vous essayez de calculer n, comme la somme suivante est égale à 1.

Comme vous pouvez le voir approche à 1 asymptotiquement, il ne sera jamais EQUAL à 1. Corrigez-moi si j'ai mal compris votre question.