Source code for factorial.pro:

;$Id: factorial.pro,v 1.12 2001/01/15 22:28:04 scottm Exp $
;
; Copyright (c) 1994-2001, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       FACTORIAL
;
; PURPOSE:
;       This function computes the factorial N! as the double-precision
;       product, (N) * (N-1) * (N-2) * ...... * 3 * 2 * 1.
;
; CATEGORY:
;       Special Functions.
;
; CALLING SEQUENCE:
;       Result = Factorial(n)
;
; INPUTS:
;       N:    A non-negative scalar or array of values.
;
; KEYWORD PARAMETERS:
;       STIRLING:    If set to a non-zero value, Stirling's asymptotic
;                    formula is used to approximate N!.
;
;       UL64: Set this keyword to return values as unsigned 64-bit integers.
;
; EXAMPLE:
;       Compute 20! with and without Stirling's asymptotic formula.
;         result_1 = factorial(20, /stirling)
;         result_2 = factorial(20)
;
;       Result_1 and result_2 should be 2.4227869e+18 and 2.4329020e+18
;       respectively.
;
; REFERENCE:
;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
;       Erwin Kreyszig
;       ISBN 0-471-55380-8
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, November 1994
;       CT, RSI, June 2000: Rewrote to handle array input; uses GAMMA for
;              non-integers; added UL64 keyword.
;       CT, RSI, October 2000: Separate calculation for scalar input.
;              Use common variable for look-up table.
;       Richard Massey, October 2003: Calls to GAMMA diverted to NR_GAMMA after
;              RSI created a built-in function GAMMA and changed its behaviour.
;-

function factorial, x, stirling = stirling, UL64=ul64

  COMPILE_OPT idl2, HIDDEN

  on_error, 2
  COMMON commonFactorial, factorialBuiltIn

	IF KEYWORD_SET(stirling) THEN BEGIN
		if MIN(x) lt 0 then $
			message, 'Values for N must be non-negative.'
		;Approximate N! using Stirling's formula.
		RETURN, SQRT(2.0d * !DPI * x) * (x / EXP(1.0d))^(x+0.0d)
	ENDIF

	IF (N_ELEMENTS(factorialBuiltIn) LT 1) THEN BEGIN
		factorialBuiltIn = $
			[1ull,1,2,6,24,120,720,5040,40320,362880,3628800, $
			39916800,479001600,6227020800,87178291200,1307674368000, $
			20922789888000,355687428096000,6402373705728000, $
			121645100408832000,2432902008176640000ull]
	ENDIF

	isInt = (x EQ FIX(x)) AND (x LE 20)
	n = N_ELEMENTS(x)

; For speed purposes, split the calculation into 1-element or array input

	IF (n EQ 1) THEN BEGIN   ; scalar or 1-element vector

		if x[0] lt 0 then $
			message, 'Values for N must be non-negative.'
		IF isInt[0] THEN BEGIN
		; do the integer factorials
			fact = factorialBuiltIn[x]
			IF NOT KEYWORD_SET(ul64) THEN fact = DOUBLE(fact)
		ENDIF ELSE BEGIN
		; do floating-point factorials using GAMMA(n+1)
		; the x*0 ensures that a 1-element vector remains a vector
			fact = x*0 + NR_GAMMA(x+1d)
			IF KEYWORD_SET(ul64) THEN fact = ULONG64(fact)
		ENDELSE

	ENDIF ELSE BEGIN   ; array

		if MIN(x) lt 0 then $
			message, 'Values for N must be non-negative.'
		fact = KEYWORD_SET(ul64) ? ULON64ARR(n) : DBLARR(n)

		; first do all the integer factorials
		whereInt = WHERE(isInt,nInt, $
			COMPLEMENT=whereNotInt, NCOMPLEMENT=nNotInt)
		IF (nInt GT 0) THEN BEGIN
			; most common cases...
			fact[whereInt] = factorialBuiltIn[x[whereInt]]
		ENDIF

		; now do all the floating-point factorials using GAMMA(n+1)
		IF (nNotInt GT 0) THEN $
			fact[whereNotInt] = NR_GAMMA(x[whereNotInt]+1d)

	ENDELSE

	RETURN, fact

end

Return to the shapelets web page or the code help menu.