Math 455.1  •  13 March 2009

The factorial function: 
An introduction to Mathematica programming

Version 7.01

Copyright © 2009 by Murray Eisenberg.  All rights reserved.

The factorial function in mathematics

Recall the meaning of the factorial function: For a positive integer n, the value n! is the product 1 · 2 · 3 · ⋯ · n of the first n positive integers, and 0! = 1. In particular, 1! = 1, 2! = 2, 3! = 6, and 4! = 24.

The purpose of this notebook is to show various ways of expressing the factorial function in Mathematica. These ways can serve as introductions to some of Mathematica's programming paradigms.

Mathematica's built-in Factorial function

The factorial function is built into Mathematica and, in its abbreviated input form, uses the very same notation as in mathematics itself:

4!

Tables of factorials

To display the values of the factorial function for a whole range of nonnegative integers, you may use the built-in function Table:

Table[n!,{n,0,10}]

For a more informative table, make each element of the resulting list not just the factorial value n!, but the list {n,n!}:

Table[{n,n!},{n,0,10}]

To display the resulting table in a nicer way, wrap the whole thing in the TableForm function (which, when used with just the Table result as argument, is more convenient to apply in the //TableForm suffix syntax:

Table[{n,n!},{n,0,10}]//TableForm

Finally, provide headings over the columns by using the TableHeadings option to TableForm:

TableForm[Table[{n,n!},{n,0,10}],TableHeadings->{None,{"n","n!"}}]

FullForm of n!

The form n! is just an abbreviation for the full form Factorial[n].  For example:

Factorial[4]

To see that that n! is in general (for nonegative integers n) an abbreviation for Factorial[n], just look at the FullForm of n!—the way Mathematica internally represents the abbreviated expression n!:

FullForm[n!]

FullForm[4!]

Now why wasn't the result here Factorial[4] rather than the computed value 24?  Because Mathematica evaluated 4! first, before even evaluating FullForm of its argument 4!. You can force Mathematica to delay evaluating the argument 4! by Holding the entire expression:

Hold[FullForm[4!]]

Factorial on lists

One of the attributes of Factorial is that it is listable …

Attributes[Factorial]

… which means that Factorial automatically distributes over an entire list of nonnegative integers given as its argument.  For example:

Factorial[{1, 3, 5}]

Or, once again using the abbreviation:

{1,3,5}!

To get a list of consecutive integers for the argument to Factorial, use the Range function:

Range[5]
Factorial[Range[5]]

Range[2,6]
Factorial[Range[2,6]]

Factorial[n] as a product

The definition n! = 1 · 2 · 3 · ⋯ · n (and 0! = 1) can be expressed by using the product notation ∏ that is analogous to the familiar sum notation ∑:
	n!=(∏_(i=1)^n)i and 0!=1.
This product notation can be entered directly into Mathematica. The easiest way to do this is by using a palette, selected from the Palettes menu. Exactly what you do depends on which version of Mathematica you are using.

Mathematica Version 7: Select the Classroom Assistant palette or the Basic Math Assistant palette. If you do not see the palette, first select it from the Palettes menu.
In the palette's Calculator section at its top, click the Advanced tab. In your own notebook, put your cursor into an input cell, such as the next input cell.
On the palette, click the ∏(Underoverscript[ , (□=□), (□   )])■ button. You'll see this in your input cell: (∏_(index=start)^end)expr
The cursor will be in the index placeholder box, which will be grayed. Now do this:
• Type i in that box.
• Press ⇥ to move the cursor to the start placeholder box and type 1 there.
• Press  ⇥ to move the cursor to the end placeholder box and type, say, 20 there.
You should now have the expression (∏_(i=1)^20)i in the cell.  Then evaluate the cell.

(* enter (∏_(i=1)^20)i in this input cell *)

Mathematica Version 6: Select the Basic Input palette. If you do not see the palette, select it from the Palettes menu. In your own notebook, put your cursor into an input cell, such as the next input cell.
On the palette, click the ∏(Underoverscript[ , (□=□), (□   )])■ button. You'll see this in your input cell: (∏_(□=□)^□)□	
The cursor will be in the leftmost little empty box—a placeholder—below the big product symbol ∏.  Now do this:
• Type i in that box.
• Press ⇥ to move the cursor to the rightmost placeholder box under the ∏ and type 1 there.
• Press  ⇥ to move the cursor to the placeholder box above the ∏ and type, say, 20 there.
You should now have the expression (∏_(i=1)^20)i in the cell.  Then evaluate the cell.

(* enter (∏_(i=1)^20)i in this input cell *)

⚠  Before proceeding, Clear any value you might have assigned to the symbol n.

Factorials by functional programming

Here's a nicer way to accomplish the same thing as using the product notation (∏_(i=1)^n)i. (It is the analog for multiplication of the  form Apply[Plus, nums] for addition that is used in Divisibility.nb.)

Apply[Times,Range[20]]

Check that the result of the preceding input cell is really the value of 20!.

You may use the @@ notation as a shortcut for using Apply:

Times@@Range[20]

Look at the FullForm of the preceding input:

FullForm[Hold[Times@@Range[20]]]

The output there confirms that func@@arg is really just a shortcut—an abbreviation—for Apply[func,arg].

Why does Apply[Times,Range[n]] have the same value as Factorial[n]?  First, what about Times?

Times[2,7]

FullForm[m n]

Hold[FullForm[2 7]]

Predict the values of each of the following:

Range[3]
Times[{1,2,3}]

Range[2]
Times[{1,2}]

Range[1]
Times[{1}]

Range[0]
Times[{}]

Let's just make sure again that the names i and n do not have any specific numeric value:

Clear[i,n]

Recall that an expression of the form fname[var_] := expr is used to define a function of your own.  For example, the expression
	f[x_]:=(x^2)
would define the squaring function.
Now define anew the factorial function in terms of Times and Apply:

factorialProduct[n_]:=Apply[Times,Range[n]]

factorialProduct[4]

Factorials by recursive programming

In this section you learn how to define the factorial function by a recursive method—a method that, in effect, expresses the factorial function in terms of itself!

The recursive definition is based upon the following two mathematical properties the factorial function:
	0!=1,                                                                           
n!=n  (n - 1)!          (n=1,2,3, …).
Apply these properties repeatedly:
	1!=1 (0!) =1 · 1=1,
2!=2 (1!)=2 · 1=2,
3!=3 (2!)=3 · 2=6,
etc.

The following definition in Mathematica is a direct realization of the two properties:

myFactorial[0]=1;
myFactorial[n_]:=n myFactorial[n-1]

It really works—try the following:

myFactorial[5]
myFactorial[20]

Another way to view what we have called a recursive definition is a definition by pattern matching. (Indeed, pattern matching is one of the underlying principles to all of Mathematica.)  The expression myFactorial[0]=1 is a particular rule. The expression myFactorial[n_]:= n myFactorial[n-1] is a general rule. So how does Mathematica know which rule to use for a given input? Mathematica always tries particular rules first. Only when no particular rule applies does it try a general rule.

The reason that those two simple rules suffice for Mathematica to calculate myFactorial[n] no matter how large n is—provided you wait long enough and your computer doesn't run out of memory—is that when Mathematica evaluates an expression, it repeatedly applies all relevant rule until it cannot further change the value.

You can follow all the steps in applying the two rules by using the Trace function.  For example:

Trace[myFactorial[4]]

If that's too much output to examine, try this version of Trace which lists only those steps in the recursion that explicitly involve the name fact but not any of the resulting arithmetic steps.

Trace[myFactorial[4],myFactorial]

The recursive definition will work with larger arguments, too. For example:

myFactorial[100]

However, Mathematica won't be able to complete the recursion if you supply an argument that is too large:

myFactorial[300]

To begin to understand why Mathematica cannot complete such a recursion, try the following and note at each evaluation what Mathematica knows about the function myFactorial.

?myFactorial

myFactorial[1]

?myFactorial

myFactorial[2]

?myFactorial

myFactorial[3]

No matter what argument you supply to myFactorial, Mathematica knows nothing more about the function than when it started. And that is why this recursive definition, as given above, is terribly inefficient. It's so inefficient that Mathematica won't even complete the calculation if you use arguments that are too large, such as in myFactorial[250].

To force Mathematica to remember values of myFactorial that it has already calculated—and thereby make the use of myFactorial much more efficient—we use the technique of caching. Here is a revised definition of myFactorial, using caching:

fact[0]=1;
fact[n_]:=fact[n]=n fact[n-1]

Try the following and note at each step what is known about the revised function fact.

?fact

fact[1]

?fact

fact[2]

?fact

fact[3]

fact[100]

?fact

fact[250]

Of course there is a downside to caching: a tradeoff of space efficiency for time efficiency. With the original, uncached version of myFactorial, a large number of calls to the function weere requred, one after the other, to find myFactorial[n] for a large n; with the revised, cached versionfact, a large number of particular rules about the function—its values at all numbers up to and including the largest used so far as argument—are accumulated and have to be stored in memory.

Factorials by recursive programming with an explicit condition

The reason that the recursive definition 
	fact[0]=1;
	fact[n_]:=n fact[n-1]
of factorial works is that Mathematica always applies particular rules before general rules, and that Mathematica repeatedly applies rules until no further changes in the result are possible.If you do not want to take full advantage of this behavior, you may use either of two alternate approaches that still provide a recursive definition—the first using the operator /; and the second using the function If.

The operator /;

To make explicit that Mathematica should use the "general" rule for factorial only when the input is positive, you may use the operator /: like this:

fac[0]:=1
fac[n_]:=n fac[n-1]/;n>0

For example:

fac[5]

You may read the rule fac[n_]:=n fac[n-1]/;n > 0 as "use the definition fac[n_]:=n fac[n-1] provided that the condition n>0 holds."

Some people like to put the test n>0 on the left-hand side of the rule, like this:

Clear[fac]
fac[0]:=1
fac[n_ /;n>0]:=n fac[n-1]

The effect here is the same as before:

fac[5]

Remember, though, that you do not actually need to use /; in this sort of simple situation, where you have a single particular rule and one general rule. (There are other, more complicated, situations, though, where /; is more useful.)

The function If

If you are fluent in some procedural programming language such as C/C++ or Java, you may feel more comfortable using Mathematica's function If to define a function whose values can take one of two different forms, depending on the form of the input.

f[n_]:=If[n==0,1,n f[n-1]]

Notice that the exression n==0  uses two equal signs to form an equation as the test (but Mathematica merges the two after you type them into what looks like a long equal sign). Think of the preceding definition as follows:
	f[n_]:=If[n==0,(* then* ) 1,(* else *) f[n-1]]
Two comments are embedded there to suggest the meaning of the second and third arguments to If:

{f[0],f[5]}

Note that, unlike in most procedural programming languages, the If in Mathematica is not a "key word" or "reserved word". Rather, If is a built-in function.

Type-checking of input

Look again at the functional programming version, factorialProduct:

?factorialProduct

What happens if you supply a type of input that was not intended to be used?

factorialProduct[-2]

factorialProduct[1/2]

factorialProduct[0.5]

What happened there, and how can it be fixed?

Every Mathematica expression has a head, which tells what kind of object it is.  For example:

{Head[-5],Head[22/3],Head[2.34],Head[β],Head[Plot[x^2,{x,0,1}]]}

At the very least, you want the input to factorialProduct to be an integer, that is, to be an expression which, when evaluated, has the Head of Integer.

{Head[5],Head[2+3],Head[-2],Head[1/2],Head[0.5]}

The way to prevent a function such as factorialProduct from evaluating its input unless that input has a Head of Integer is to append the desired Head to the input variable pattern n_, like this:

factorialProduct[n_Integer]:=Apply[Times,Range[n]]

factorialProduct[5]
factorialProduct[1/2]
factorialProduct[0.5]

Oops! Why didn't the new definition do what was expected now? Look at the current rules associated with the name factorialProduct:

?factorialProduct

There are two rules: the more particular one just set down that restricts the input to be an integer; and the more general one from earlier that accepts any input! Mathematica always tries the more particular rules first and, if it cannot apply those to evaluate the given expression, then it tries the more general rules. So either you have to use a new name for the intended new function or else Clear the rules associated with the old function.  Let's do the latter:

Clear[factorialProduct]
factorialProduct[n_Integer]:=Apply[Times,Range[n]]

factorialProduct[5]
factorialProduct[1/2]
factorialProduct[0.5]

So far OK. But what about the following now?

factorialProduct[-2]

A new way is needed to prevent factorialProduct[n] from being evaluated when n is not a nonnegative integer.  Just examining the Head of the input will not suffice:

Head[-2]

Instead, a new kind of restriction is needed on the input pattern n_. One way to impose such a restriction is to append an expression of the form ?test, where test is a function that takes only the values True and False. Here we use the function NonNegative:

Clear[factorialProduct]
factorialProduct[n_?NonNegative]:=Apply[Times,Range[n]]

factorialProduct[5]
factorialProduct[-2]

factorialProduct[1/2]

Oops! The lastest version of factorialProduct checks whether the input is nonnegative, but it no longer checks whether it is integer. We need to check for both types:

Clear[factorialProduct]
factorialProduct[n_Integer?NonNegative]:=Apply[Times,Range[n]]

factorialProduct[5]
factorialProduct[-2]
factorialProduct[1/2]

All is well!

Appendix: Factorials by iterative programming

In this section you learn how to define the factorial function by an iterative procedure. Of course such an iterative procedure superfluous—because Factorial is a built-in function, and because even if Factorial were not built-in, you could define it by means of the simple "functional programming" expression involving Apply and Times.

First use a Do expression.

?Do

Notice that Do is a full-fledged function, not just a "keyword". However, it returns a Null result that ordinarily does not display in an Out cell. To see the final result found by a Do expression, you have to print or otherwise evaluate that final value.

Do[(i^2),{i, 5}]

Do[Print[(i^2)], {i,5}]

Note that the iterator {i,5} is shorhand for {i,1,5}.

Do[temp=(i^2),{i,5}];temp

OK, now use Do to calculate Factorial[5]:

runningProduct = 1;
Do[runningProduct=i×runningProduct,{i,5}];
runningProduct

The three separate expressions in the preceding input cell, separated by semicolons, constitute a compound expression.  To make such a compound expression calculate Factorial[n] for arbitrary n and not just the particular n=5, change 5 to n and encapsulate the entire compound expression in a Module, like this:

factorialDo[n_]:=Module[{runningProduct,i},   (* declare local variables *)
runningProduct =1;
Do[runningProduct=i×runningProduct,{i,n}];
runningProduct
]

If you wish, you may initialize the local variable runningProduct within the local variable declaration:

Clear[factorialDo]
factorialDo[n_]:=Module[{runningProduct=1,i},
Do[runningProduct=i×runningProduct,{i,n}];
runningProduct
]

factorialDo[500]

If it makes you more comfortable, you can include in the procedural program factorialDo a separate test for input 0—even though, as you should have observed, such a test is really superfluous.  Here is the (unnecessarily) modified form:

overlyCautiousFactorialDo[n_] :=  Module[{fact,i},
If[n == 0, (* then *) 1, 
(* else *)
	fact = 1;
	Do[fact = i × fact, {i, n}];
fact]
]

overlyCautiousFactorialDo[0]