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.

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>

Numerical Approximation

<Ada> with Ada.Numerics.Generic_Complex_Types; with Ada.Numerics.Generic_Complex_Elementary_Functions; with Ada.Numerics.Generic_Elementary_Functions; with Ada.Text_IO.Complex_Io; with Ada.Text_Io; use Ada.Text_Io;

procedure Factorial_Numeric_Approximation is

  type Real is digits 15;
  package Complex_Pck is new Ada.Numerics.Generic_Complex_Types(Real);
  use Complex_Pck;
  package Complex_Io is new Ada.Text_Io.Complex_Io(Complex_Pck);
  use Complex_IO;
  package Cmplx_Elem_Funcs is new Ada.Numerics.Generic_Complex_Elementary_Functions(Complex_Pck);
  use Cmplx_Elem_Funcs;
  
  function Gamma(X : Complex) return Complex is
     package Elem_Funcs is new Ada.Numerics.Generic_Elementary_Functions(Real);
     use Elem_Funcs;
     use Ada.Numerics;
     -- Coefficients used by the GNU Scientific Library
     G : Natural := 7;
     P : constant array (Natural range 0..G + 1) of Real := (
        0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7);
     Z : Complex := X;
     Cx : Complex;
     Ct : Complex;
  begin
     if Re(Z) < 0.5 then
        return Pi / (Sin(Pi * Z) * Gamma(1.0 - Z));
     else
        Z := Z - 1.0;
        Set_Re(Cx, P(0));
        Set_Im(Cx, 0.0);
        for I in 1..P'Last loop
           Cx := Cx + (P(I) / (Z + Real(I)));
        end loop;
        Ct := Z + Real(G) + 0.5;
        return Sqrt(2.0 * Pi) * Ct**(Z + 0.5) * Exp(-Ct) * Cx;
     end if;
  end Gamma;
  
  function Factorial(N : Complex) return Complex is
  begin
     return Gamma(N + 1.0);
  end Factorial;
  Arg : Complex;

begin

  Put("factorial(-0.5)**2.0 = ");
  Set_Re(Arg, -0.5);
  Set_Im(Arg, 0.0);
  Put(Item => Factorial(Arg) **2.0, Fore => 1, Aft => 8, Exp => 0);
  New_Line;
  for I in 0..9 loop
     Set_Re(Arg, Real(I));
     Set_Im(Arg, 0.0);
     Put("factorial(" & Integer'Image(I) & ") = ");
     Put(Item => Factorial(Arg), Fore => 6, Aft => 8, Exp => 0);
     New_Line;
  end loop;

end Factorial_Numeric_Approximation; </Ada> Output:

factorial(-0.5)**2.0 = (3.14159265,0.00000000)
factorial( 0) = (     1.00000000,     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)

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
;

bash

fact ()
{
  if [ $1 -eq 0 ];
    then echo 1;
    else echo $(($1 * $(fact $(($1-1)) ) ));
  fi;
}

C

Iterative

<c> int factorial(int n) {

 int result = 1;
 for (int i = 1; i <= n; ++i)
   result *= i;
 return result;

} </c>

Recursive

<c> int factorial(int n) {

 if (n == 0)
   return 1;
 else
   return n*factorial(n-1);

} </c>

C++

The C versions work unchanged with C++, however, here is another possibility using the STL and boost: <cpp>

1. include <boost/config.hpp>
2. include <boost/counting_iterator.hpp>

int factorial(int n) {

 boost::counting_iterator_generator<int>::type first(1), last(n+1); // last is one-past-end
 return std::accumulate(first, last, 1, std::multiplies<int>());

} </cpp>

Forth

: fac ( n -- n! ) 1 swap 1+ 1 ?do i * loop ;

Fortran

Works with: Fortran version 90 and later

A simple one-liner is sufficient

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

Works with: Fortran version 77 and later

      INTEGER FUNCTION FACT(N)
      INTEGER N, I
      FACT = 1
      DO 10, I = 1, N
        FACT = FACT * I
 10   CONTINUE
      RETURN
      END

Haskell

The simplest description: factorial is the product of the numbers from 1 to n:

factorial n = product [1..n]

Or, written explicitly as a fold:

factorial n = foldl (*) 1 [1..n]

See also: The Evolution of a Haskell Programmer

Java

Iterative

<java>public static long fact(int n){

  if(n < 0){
     System.err.println("No negative numbers");
     return 0;
  }
  long ans = 1;
  for(int i = 1;i <= n;i++){
     ans *= i;
  }
  return ans;

}</java>

Recursive

<java>public static long fact(int n){

  if(n < 0){
     System.err.println("No negative numbers");
     return 0;
  }
  if(n == 0) return 1;
  return n * fact(n-1);

}</java>

JavaScript

<javascript>function factorial(n) {

 return n<2 ? 1 : n * factorial(n-1);

}</javascript>

Logo

to factorial :n
  if :n < 2 [output 1]
  output :n * factorial :n-1
end

Mathematica

Note that Mathematica already comes with a factorial function, which can be used as e.g. 5! (gives 120). So the following implementations are only of pedagogical value.

Recursive

factorial[n_Integer] := n*factorial[n-1]
factorial[0] = 1

Iterative (direct loop)

factorial[n_Integer] := 
  Block[{i, result = 1}, For[i = 1, i <= n, ++i, result *= i]; result]

Iterative (list)

factorial[n_Integer] := Block[{i}, Times @@ Table[i, {i, n}]]

OCaml

Recursive

<ocaml>let rec factorial n =

 if n <= 0 then 1
 else n * factorial (n-1)</ocaml>

The following is tail-recursive, so it is effectively iterative: <ocaml>let factorial n =

 let rec loop i accum =
   if i > n then accum
   else loop (i + 1) (accum * i)
 in loop 1 1</ocaml>

Iterative

It can be done using explicit state, but this is usually discouraged in a functional language: <ocaml>let factorial n =

 let result = ref 1 in
 for i = 1 to n do
   result := !result * i
 done;
 !result</ocaml>

Pascal

Iterative

<pascal> function factorial(n: integer): integer;

var
 i, result: integer;
begin
 result := 1;
 for i := 1 to n do
  result := result * i;
 factorial := result
end;

</pascal>

Recursive

<pascal> function factorial(n: integer): integer;

begin
 if n = 0
  then
   factorial := 1
  else
   factorial := n*factorial(n-1)
end;

</pascal>

Perl

Iterative

<perl> sub factorial {

 my $n = shift;
 my $result = 1;
 for (my $i = 1; $i <= $n; ++$i)
 {
   $result *= $i;
 };
 $result;

} </perl>

Recursive

<perl> sub factorial {

 my $n = shift;
 ($n == 0)? 1 : $n*factorial($n-1);

} </perl>

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>

Scheme

Recursive

<scheme>(define (factorial n)

 (if (<= n 0)
     1
     (* n factorial (- n 1))))</scheme>

The following is tail-recursive, so it is effectively iterative: <scheme>(define (factorial n)

 (let loop ((i 1)
            (accum 1))
   (if (> i n)
       accum
       (loop (+ i 1) (* accum i)))))</scheme>

Iterative

<scheme>(define (factorial n)

 (do ((i 1 (+ i 1))
      (accum 1 (* accum i)))
     ((> i n) accum)))</scheme>