Factorial

Factorial
You are encouraged to solve this task according to the task description, using any language you may know.

The Factorial Function of a positive integer, n, is defined as the product of the sequence n, n-1, n-2, ...1 and the factorial of zero, 0, is [defined] as being 1.

Write a function to return the factorial of a number. Solutions can be iterative or recursive (though recursive solutions are generally considered too slow and are mostly used as an exercise in recursion). Support for trapping negative n errors is optional.

References

Wikipedia

Ada

Iterative

<Ada> function Factorial (N : Positive) return Positive is

  Result : Positive := N;
  Counter : Natural := N - 1;

begin

  for I in reverse 1..Counter loop
     Result := Result * I;
  end loop;
  return Result;

end Factorial; </Ada>

Recursive

<Ada> function Factorial(N : Positive) return Positive is

  Result : Positive := 1;

begin

  if N > 1 then
     Result := N * Factorial(N - 1);
  end if;
  return Result;

end Factorial; </Ada>

ALGOL 68

Iterative

PROC factorial = (INT upb n)LONG LONG INT:(
  LONG LONG INT z := 1;
  FOR n TO upb n DO z *:= n OD;
  z
);

Numerical Approximation

# Coefficients used by the GNU Scientific Library #
INT g = 7;
[]REAL p = []REAL(0.99999999999980993, 676.5203681218851, -1259.1392167224028, 
  771.32342877765313, -176.61502916214059, 12.507343278686905, 
  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)[@0];

PROC complex gamma = (COMPL in z)COMPL: (
  # Reflection formula #
  COMPL z := in z;
  IF re OF z < 0.5 THEN
    pi / (complex sin(pi*z)*complex gamma(1-z))
  ELSE
    z -:= 1;
    COMPL x := p[0];
    FOR i TO g+1 DO x +:= p[i]/(z+i) OD;
    COMPL t := z + g + 0.5;
    complex sqrt(2*pi) * t**(z+0.5) * complex exp(-t) * x
  FI
);

OP ** = (COMPL z, p)COMPL: ( z=0|0|complex exp(complex ln(z)*p) );
PROC factorial = (COMPL n)COMPL: complex gamma(n+1);

FORMAT compl fmt = $g(-16, 8)"⊥"g(-10, 8)$;
printf(($q"factorial(-0.5)**2="f(compl fmt)l$, factorial(-0.5)**2));
FOR i TO 9 DO
  printf(($q"factorial("d")="f(compl fmt)l$, i, factorial(i)))
OD

Output:

factorial(-0.5)**2=      3.14159265⊥0.00000000
factorial(1)=      1.00000000⊥0.00000000
factorial(2)=      2.00000000⊥0.00000000
factorial(3)=      6.00000000⊥0.00000000
factorial(4)=     24.00000000⊥0.00000000
factorial(5)=    120.00000000⊥0.00000000
factorial(6)=    720.00000000⊥0.00000000
factorial(7)=   5040.00000000⊥0.00000000
factorial(8)=  40320.00000000⊥0.00000000
factorial(9)= 362880.00000000⊥0.00000000

Recursive

PROC factorial = (INT n)LONG LONG INT:
  CASE n+1 IN
    1,1,2,6,24,120,720 # a brief lookup #
  OUT
    n*factorial(n-1)
  ESAC
;

Fortran

Works with: Fortran version 90 and later

A simple one-liner is sufficient

nfactorial = PRODUCT((/(i, i=1,n)/))

Python

Iterative

<python>def factorial(n):

   if n == 0:
       return 1
   z=n
   while n>1:
       n=n-1
       z=z*n
   return z

</python>

Functional

<python>from operator import mul

def factorial(n):

   return reduce(mul, xrange(1,n+1), 1)

</python>

Sample output:

<python>

>>> for i in range(6):

print i, factorial(i)

0 1

1 1

2 2

3 6

4 24

5 120

>>>

</python>

Numerical Approximation

The following sample uses Lanczos approximation from http://en.wikipedia.org/wiki/Lanczos_approximation <python> from cmath import *

1. Coefficients used by the GNU Scientific Library

g = 7 p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,

    771.32342877765313, -176.61502916214059, 12.507343278686905,
    -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]

def gamma(z):

 z = complex(z)
 # Reflection formula
 if z.real < 0.5:
   return pi / (sin(pi*z)*gamma(1-z))
 else:
   z -= 1
   x = p[0]
   for i in range(1, g+2):
     x += p[i]/(z+i)
   t = z + g + 0.5
   return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x

def factorial(n):

 return gamma(n+1)

print "factorial(-0.5)**2=",factorial(-0.5)**2 for i in range(10):

 print "factorial(%d)=%s"%(i,factorial(i))

</python> Output:

factorial(-0.5)**2= (3.14159265359+0j)
factorial(0)=(1+0j)
factorial(1)=(1+0j)
factorial(2)=(2+0j)
factorial(3)=(6+0j)
factorial(4)=(24+0j)
factorial(5)=(120+0j)
factorial(6)=(720+0j)
factorial(7)=(5040+0j)
factorial(8)=(40320+0j)
factorial(9)=(362880+0j)

Recursive

<python>def factorial(n):

   z=1
   if n>1:
       z=n*factorial(n-1)
   return z

</python>