Factoriel et log factoriel

Factoriel:

Quand vous avez besoin de calculer n!, il existe plusieurs solutions: 

1. La solution « rapide »: qui utilise une boucle ou une fonction récursive: 
   

   def factorial_for(n): r = 1 for i in range(2, n + 1): r *= i return(r)

   def factorial_rec(n): if n > 1: return(n * factorial_rec(n - 1)) else: return(1)

   Ici, la multiplication séquentielle de chaque nombre va générer un nombre très grand très rapidement, ce qui n’est pas optimum. En effet, les ordinateurs sont plus rapides quand ils multiplient des petits nombres (120×30240) par rapport à la même multiplication avec un grand nombre (362880×10).

2. La solution « intelligente »: Ici, l’idée est de séparer une séquence de nombres à multiplier en deux par le milieu, de façon récursive.
   

   10!	= 1x2x3x4x5x6x7x8x9 x 10	= 362880 x 10	(rush solution)
   	= 1x2x3x4x5 x 6x7x8x9x10	= 120 x 30240	(1st split)
   	= 1x2x3 x 4×5 x 6x7x8 x 9×10	= 6 x 20 x 336 x 90	(2nd split)

   def range_prod(lo, hi):
     if lo + 1 < hi:
       mid = (hi + lo) // 2
       return(range_prod(lo, mid) * range_prod(mid + 1, hi))
     if lo == hi:
       return(lo)
     return(lo * hi)
   

   def factorial_tree(n):
     if n < 2:
       return(1)
     return(range_prod(1, n))

3. La solution « paresseuse » ou « très intelligente » (tout dépend du point de vue): consiste à calculer une estimation du factoriel avec l’approximation de Stirling. 

   import math
   def stirling(n):
     return(math.sqrt(2 * math.pi * n) * (n / math.e) ** n)

   Cette implémentation est limitée à n<170. Toutefois, comme vous pouvez le constater plus bas, avec n=170, le résultat est énorme.
   Avec ce genre de gros nombre, nous préférons de toute façon travailler avec les log factoriels.

   In [01]: stirling(170)
   Out[01]: 7.253858934542908e+306
   In [02]: stirling(171)
   Out[02]: inf
   

    

Log Factorial:

Pour les deux premières solutions, vous n’avez qu’à ajouter math.log() (en faisant import math d’abord) :

def log_fact_rec(n):
  return(math.log(factorial_rec(n)))

def log_fact_tree(n):
  return(mat.log(factorial_tree(n)))

def log_fact_for(n):
  return(math.log(factorial_for(n)))

Le log factoriel avec l’approximation de Stirling est implémentée comme ceci :

def log_stirling(n):
 fn = float(n)
 fn = n * math.log(fn) + math.log(2 * math.pi * fn) / 2 - fn + \
      (fn ** -1) / 12 - (fn ** -3) / 360 + (fn ** -5) / 1260 - \
      (fn ** -7) / 1680 + (fn ** -9) / 1188
 return(fn)

Note:

Dans le cas où vous avez de multiples log factoriels à calculer, il peut être utile de sauvegarder et réutiliser des résultats précédents et ainsi gagner du temps (par exemple, 10! = 6!x7x8x9x10) comme ceci :

def log_fact_exact_saved(n, value_saved):
  select = max((i for i in value_saved.keys() if i <= n))
  if select != n:
    value_saved[n] = range_prod(select + 1, n) * value_saved[select]
  return(math.log(value_saved[n]))

def log_fact_approx_saved(n, value_saved):
  if n not in value_saved:
    value_saved[n] = log_stirling(n)
  return(value_saved[n])

Test:

Nous pouvons maintenant tester et comparer chacune de ces solutions en utilisant les librairies math et scipy pour voir laquelle est la plus rapide :

import sys
sys.setrecursionlimit(1000000) // increase the depth 
import math
import scipy.special
def test_saved(saved, approx):
    value = [10, 50, 100, 1000, 1000, 2000, 5000, 5000, 5000,
             5000, 10000, 11000, 11000, 11000, 20000, 20000]
    if approx:
        if saved:
            value_saved = {1: 1}
            for num in value:
                log_fact_approx_saved(num, value_saved)
        else:
            for num in value:
                log_stirling(num)
    else:
        if saved:
            value_saved = {1: 1}
            for num in value:
                log_fact_exact_saved(num, value_saved)
        else:
            for num in value:
                log_fact_tree(num)

def log_fact_math(n):
    return(math.log(math.factorial(n)))

def log_fact_scipy(n, exact=False):
    return(math.log(scipy.special.factorial(n, exact)))

import timeit
num = 10000
rep = 200
print("rec:" + str(timeit.timeit("log_fact_rec(" + str(num) + ")",
      setup="from __main__ import log_fact_rec", number=rep)))
print("for:" + str(timeit.timeit("log_fact_for(" + str(num) + ")",
      setup="from __main__ import log_fact_for", number=rep)))
print("scipy exact:" + str(timeit.timeit("log_fact_scipy(" + str(num) + ", True)",
      setup="from __main__ import log_fact_scipy", number=rep)))
print("math:" + str(timeit.timeit("log_fact_math(" + str(num) + ")",
      setup="from __main__ import log_fact_math", number=rep)))
print("tree:" + str(timeit.timeit("log_fact_tree(" + str(num) + ")",
      setup="from __main__ import log_fact_tree", number=rep)))

print("scipy approx:" + str(timeit.timeit("log_fact_scipy(" + str(num) + ")",
      setup="from __main__ import log_fact_scipy", number=rep)))
print("log_stirling:" + str(timeit.timeit("log_stirling(" + str(num) + ")",
      setup="from __main__ import log_stirling", number=rep)))

print("value log(100!) tree:" + str(log_fact_tree(100)))
print("value log(100!) scipy approx:" + str(log_fact_scipy(100)))
print("value log(100!) log_stirling:" + str(log_stirling(100)))

print("severals factorial:" + str(timeit.timeit("test_saved(False, False)",
      setup="from __main__ import test_saved", number=rep)))
print("severals factorial saved:" + str(timeit.timeit("test_saved(True, False)",
      setup="from __main__ import test_saved", number=rep)))
print("severals factorial approx:" + str(timeit.timeit("test_saved(False, True)",
      setup="from __main__ import test_saved", number=rep)))
print("severals factorial approx saved:" + str(timeit.timeit("test_saved(True, True)",
      setup="from __main__ import test_saved", number=rep)))

En utilisant python 2.7.6 et scipy 0.17.0, nous avons :

J’espère que cet article vous aidera à choisir la meilleure solution (adaptée à vos besoins!) pour calculer des nombres factoriels et log factoriels.