Chapter 1. Numerics

Jenny I’ve got your number I need to make you mine Jenny don’t change your number Eight six seven five three oh nine Eight six seven five three oh nine Eight six seven five three oh nine Eight six seven five three oh nine

Mathematics is a huge, almost all-encompassing subject, and the average layperson often fails to appreciate the types of exotic objects that are in the mathematician’s domain. Yet every person on the street perceives math is about numbers. So even though numbers only scratch the surface of math and Mathematica, it makes sense to begin with their representation.

Mathematica supports four numerical types: Integer, Rational, Real, and Complex. In the following examples we use Mathematica’s comment notation ( *comment* ).

1         (*The integer one*)
1 / 2     (*The rational one half*)
1.2 * ^ 8 (*The real 1.2 x 10^8*)
3 + 2 I   (*The complex number 3+2i*)

There is no need to take my word that these expressions have the specified types. You can ask Mathematica to tell you using the function Head[] , which returns the head of an expression (i.e., head of a list).

In[2]:=  Head[1]
Out[2]=  Integer

In[3]:=  Head[1/2]
Out[3]=  Rational

In[4]:=  Head[1.2 ^ 8]
Out[4]=  Real

In[5]:=  Head[3 + 2 I]
Out[5]=  Complex

Although Mathematica does not internally store numbers as lists, it provides the illusion that a number has a head indicating its type. This is consistent with the fact that everything in Mathematica is an expression and every expression must have a head. It is also common for Mathematica to use the head to indicate type when constructing more complex objects. See 1.5 Working with Intervals, for example. If you are confused by this, for now, just think of Head as returning a type name when presented with an atomic expression (expressions that can’t be divided into subexpressions).

Mathematica is unique in comparison to most mathematical tools and programming languages in that it will usually produce exact results unless you tell it otherwise. The following examples show the difference between exact and approximate results. 1.1 Controlling Precision and Accuracy and 1.2 Mixing Different Numerical Types show you how to make Mathematica use the appropriate form.

Exact results are displayed in their entirety when possible or symbolically when full display would be impossible due to the infinity of the exact representation.

Exact and Approximate Results

Approximate numeric results are represented in machine precision floating point by default. On most modern computers, this means double-precision floating-point numbers, which contain a total of 64 binary bits, typically yielding 16 decimal digits of mantissa. You can also specify numbers with greater than machine precision (see 1.1 Controlling Precision and Accuracy) but there is a performance cost: Mathematica must switch from the native hardware-based floating-point algorithms to software-based ones.

In[8]:= 3. ^ 1000
Out[8]= 1.322070819480807 x 10477

In[9]:= Sqrt[2.]
Out[9]= 1.41421

By adding a decimal point to a number, you force Mathematica to treat it as approximate. These approximate numbers will be machine precision by default, but there are several ways to force higher precision. 1.1 Controlling Precision and Accuracy and 1.2 Mixing Different Numerical Types in this chapter will elaborate on these differences.

For most purposes, treat precision as the total number of digits in the decimal representation of a number and accuracy as the total number of digits after the decimal. As such, precision is a measure of relative uncertainty (given a precision p a larger number will have more uncertainty than a smaller number). Accuracy is an absolute measure of uncertainty because the number of places after the decimal is independent of the magnitude of the number. Typically you only need to control precision in most applications.

There are two common syntaxes for using N[]. You already saw the functional syntax in the solution section. The second uses Mathematica’s postfix notation. See the sidebar Mathematica Expressions for a discussion of postfix and other notations.

In[29]:= Sqrt[2] //N
Out[29]= 1.41421

It is common to use this notation to force Mathematica to convert an exact or symbolic result to an approximate result as the last step in a computation. When you use postfix notation, you can explicitly specify the precision, but it is a bit awkward.

In[30]:= Sqrt[2] //N[#, 10]&
Out[30]= 1.414213562

When you don’t specify precision, Mathematica uses MachinePrecision, which is a built-in symbol that denotes the precision native to your computer’s floating-point capabilities. The numerical value of MachinePrecision is stored in a variable $MachinePrecision.

In[31]:= $MachinePrecision
Out[31]= 15.9546

There is another notation that is less common but you may come across it in Mathematica output. If a literal number is displayed with a trailing ` (backtick) followed optionally by a number, this indicates the number is either in machine precision or is in the precision specified by the number following the backtick.

In[32]:= 20` (*20 in machine precision*)
Out[32]= 20.

In[33]:= 20`20 (*20 with high precision of 20 digits*)
Out[33]= 20.000000000000000000

In a complex expression with a lot of high-precision numbers, you can avoid specifying each precision individually by using SetPrecision[].

In[34]:= SetPrecision[20. + 1/3 * 12.3 / 37.8 + Pi, 20]
         (*All numbers will be set to a precision of 20.*)
Out[34]= 23.250058262055400604

If you have an expression and need to know the precision or accuracy, you can use the following functions.

In[35]:= Precision[2.]
Out[35]= MachinePrecision

In[36]:= Precision[2'20]
Out[36]= 20.

Exact results have infinite precision.

In[37]:= Precision[Sqrt[2]]
Out[37]= ∞

In[38]:= Precision[Sqrt[2.]]
Out[38]= MachinePrecision

In[39]:= Accuracy[2.]
Out[39]= 15.6536

You are not guaranteed the accuracy you specify if the precision is too small.

In[40]:= Accuracy[N[30, {20, 20}]]
Out[40]= 18.5229

With enough precision, however, you will get accuracy.

In[41]:= Accuracy[N[30, {30, 20}]]
Out[41]= 20.

And precision can even be specified as infinite!

In[42]:= Accuracy[N[30, {Infinity, 20}]]
Out[42]= 20.

Mathematica also defines two internal variables: $MinPrecision, whose default value is zero, and $MaxPrecision, whose default value is plus infinity.

In[43]:= {$MinPrecision, $MaxPrecision}
Out[43]= {0, ∞}

You can control precision within a complex calculation (without using N[] on every intermediate result) by changing these values; however, you should only do so within a Block (a local context). For example, compare the difference between a calculation with automatic precision for intermediate results to the same calculation with fixed precision (obtained by making $MinPrecision == $MaxPrecision). Note that we must still start out the calculation with base values of at least $MinPrecision, otherwise the value will revert to the lowest precision, as explained in 1.2 Mixing Different Numerical Types.

In[44]:=  SetPrecision[(1 + Exp[Sqrt[2] + Sqrt[3]]) / 2^25, 32]
Out[44]=  7.226780742612584668840452114476x10–7

In[45]:=  Block[{$MinPrecision = 32, $MaxPrecision = 32},
           SetPrecision[(1 + Exp[Sqrt[2] + Sqrt[3]])/2^25, 32]]
Out[45]=  7.2267807426125846688404521144759x10–7

However, unless you have a very specific reason to control precision yourself, it is generally best to let Mathematica automatically handle this for you.

The Wolfram documentation for N[] is here: http://bit.ly/XVe2E.

Discussions of precision and accuracy can be found at http://bit.ly/15qq2N and http://bit.ly/icrh1 .

The most thorough discussions of precision and accuracy in Mathematica can be found in Chapter 8 of An Introduction to Programming with Mathematica (Cambridge University Press) and The Mathematica GuideBook for Numerics (Springer).

A nice essay by David Goldberg called "What Every Computer Scientist Should Know About Floating-Point Arithmetic" can be found at http://bit.ly/1EJ23y .

The most thorough discussions of Mathematica’s numerical rules can be found in Chapter 8 of An Introduction to Programming with Mathematica and The Mathematica GuideBook for Numerics.

1.5 Working with Intervals shows how to extract digits of a number in alternate bases.

Both RealDigits[] and IntegerDigits[] take the desired base and the number of desired digits (length) as optional second and third arguments, respectively.

   In[69]:=  12 !
   Out[69]=  479 001 600

   In[70]:=  IntegerDigits[12!, 10, 5]
   Out[70]=  {0, 1, 6, 0, 0}

   In[71]:=  12! // BaseForm[#, 16] & (*Consider 12! in base 16.*)
Out[71]//BaseForm=
             1c8cfc0016

   In[72]:=  IntegerDigits[12!, 16]   (*Notice how IntegerDigits
              with base 16 gives the digit values in base 10.*)
   Out[72]=  {1, 12, 8, 12, 15, 12, 0, 0}

   In[73]:=  IntegerDigits[12!, 16] // BaseForm[#, 16]&
              (*But you can easily force them to base 16.*)
Out[73]//BaseForm=
             {116, c16, 816, c16, f16, c16, 016, 016}

RealDigits can take an additional fourth argument that specifies where in the decimal expansion to start. If b is the base, then the fourth argument n means to start the counting at the coefficient signified by b^n . The following examples should clarify.

In[74]:=  N[Pi, 10]      (*Pi to 10 digits of precision.*)
Out[74]=  3.141592654

In[75]:=  RealDigits[Pi, 10, 3]
          (*Extract first three digits. Decimal place is indicated as 1.*)
Out[75]=  {{3, 1, 4}, 1}

Start at 10^-2 = 0.01, or the second digit after the decimal.

In[76]:=  RealDigits[Pi, 10, 3, -2]
          (*Extract third to fifth digit. Decimal place is indicated as -2.*)
Out[76]=  {{4, 1, 5}, -1}

Start at 10^-5 = 0.00001, or the fifth digit after the decimal.

   In[77]:=  RealDigits [Pi, 10, 3, -5]
   Out[77]=  {{9, 2, 6}, -4}

   In[78]:=  N[Pi,10] // BaseForm[#, 2] &
Out[78]//BaseForm=
             11.00100100001111110110101010001002

Here we get the digits of pi in base 2.

In[79]:=  RealDigits[Pi, 2, 5, -2]
Out[79]=  {{0, 1, 0, 0, 1}, -1}

Here is an interesting application in which IntegerDigits is combined with the Tuples function and a bit of pattern matching to get all n digits without calling IntegerDigits[] more than once. We used Short to elide the full list. (Short places <<n>> in the output to indicate n missing items.)

   In[80]:=  Tuples[IntegerDigits[43210], 4] // Short[#, 4] &
Out[80]//Short=
             {{4, 4, 4, 4}, {4, 4, 4, 3}, {4, 4, 4, 2}, {4, 4, 4, 1}, {4, 4, 4, 0}, {4, 4, 3, 4},
              {4, 4, 3, 3}, {4, 4, 3, 2}, {4, 4, 3, 1}, {4, 4, 3, 0}, {4, 4, 2, 4},
              {4, 4, 2, 3}, {4, 4, 2, 2}, {4, 4, 2, 1}, {4, 4, 2, 0}, {4, 4, 1, 4},
              {4, 4, 1, 3}, {4, 4, 1, 2}, {4, 4, 1, 1}, {4, 4, 1, 0}, {4, 4, 0, 4},
              {4, 4, 0, 3}, {4, 4, 0, 2}, <<579>>,   {0, 0, 4, 2}, {0, 0, 4, 1}, {0, 0, 4, 0},
              {0, 0, 3, 4}, {0, 0, 3, 3}, {0, 0, 3, 2}, {0, 0, 3, 1}, {0, 0, 3, 0},
              {0, 0, 2, 4}, {0, 0, 2, 3}, {0, 0, 2, 2}, {0, 0, 2, 1}, {0, 0, 2, 0},
              {0, 0, 1, 4}, {0, 0, 1, 3}, {0, 0, 1, 2}, {0, 0, 1, 1}, {0, 0, 1, 0},
              {0, 0, 0, 4}, {0, 0, 0, 3}, {0, 0, 0, 2}, {0, 0, 0, 1}, {0, 0, 0, 0}}

If you do not want the cases with leading zeros, you can use DeleteCases as follows.

   In[81]:=  DeleteCases[Tuples[IntegerDigits[43210], 4],
               {z__ /; z == 0, n__}] //  Short[#, 4] &
Out[81]//Short=
             {{4, 4, 4, 4}, {4, 4, 4, 3}, {4, 4, 4, 2}, {4, 4, 4, 1}, {4, 4, 4, 0}, {4, 4, 3, 4},
              {4, 4, 3, 3}, {4, 4, 3, 2}, {4, 4, 3, 1}, {4, 4, 3, 0}, {4, 4, 2, 4},
              {4, 4, 2, 3}, {4, 4, 2, 2}, {4, 4, 2, 1}, {4, 4, 2, 0}, {4, 4, 1, 4},
              {4, 4, 1, 3}, {4, 4, 1, 2}, {4, 4, 1, 1}, {4, 4, 1, 0}, {4, 4, 0, 4},
              {4, 4, 0, 3}, {4, 4, 0, 2}, <<454>>, {1, 0, 4, 2}, {1, 0, 4, 1}, {1, 0, 4, 0},
              {1, 0, 3, 4}, {1, 0, 3, 3}, {1, 0, 3, 2}, {1, 0, 3, 1}, {1, 0, 3, 0},
              {1, 0, 2, 4}, {1, 0, 2, 3}, {1, 0, 2, 2}, {1, 0, 2, 1}, {1, 0, 2, 0},
              {1, 0, 1, 4}, {1, 0, 1, 3}, {1, 0, 1, 2}, {1, 0, 1, 1}, {1, 0, 1, 0},
              {1, 0, 0, 4}, {1, 0, 0, 3}, {1, 0, 0, 2}, {1, 0, 0, 1}, {1, 0, 0, 0}}

The inverse of IntegerDigits[] is FromDigits[].

In[82]:= FromDigits[IntegerDigits[987654321]]
Out[82]= 987 654 321

In[83]:= FromDigits[IntegerDigits[987654321, 2], 2] (*Base 2*)
Out[83]= 987 654 321

FromDigits[] has the added capability of converting strings and roman numerals.

In[84]:= FromDigits["4750"] + 1
Out[84]= 4751

In[85]:= FromDigits["MMXIX", "Roman"] – 10
Out[85]= 2009

IntegerString[] is used to convert back to string form. I use InputForm only so the quotes are displayed.

   In[86]:=  IntegerString[4750] //InputForm
Out[86]//InputForm=
             "4750"

   In[87]:=  IntegerString[2009, "Roman"] // InputForm
Out[87]//InputForm=
             "MMIX"

On the surface, the solutions here are rather simple. In day-to-day usage, numeric conversion will not present many challenges. However, there are subtle issues and interesting theory underlying the apparent simplicity. Let’s consider rounding. Suppose you need to round a set of numbers, but the numbers still must satisfy some constraint after the rounding. Consider percentages or probabilities. One would want percentages to still add to 100 and probabilities to still sum to 1. Another context is in statistics, where we want to round while preserving certain statistical properties, such as the variance. Various forms of stochastic rounding can be used in these cases. One form of stochastic rounding that gives good results is the unbiased rounding rule. According to this rule, a number of the form x.v is rounded up with the probability v /10 and rounded down with probability (10-v)/10. So, for example, 10.5 would have equal probability of going to 10 as to 11, whereas 10.85 would have probability of 0.85 of rounding up and 0.15 of rounding down.

In[120]:= UnbiasedRound[x_] := Block[{whole = Floor[x], v},
            v = 10 * (x - whole); whole + Floor[v/10 + RandomReal[]]]

In[121]:= Table[UnbiasedRound[10.5], {20}]
Out[121]= {11, 11, 10, 11, 10, 10, 10, 11, 11, 11, 10, 11, 11, 10, 10, 11, 11, 11, 11, 11}

In[122]:= Table[UnbiasedRound[10.1], {20}]
Out[122]= {10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 10, 10, 10, 10, 11, 10, 10, 10}

In[123]:= Table[UnbiasedRound[10.8], {20}]
Out[123]= {11, 11, 11, 10, 11, 11, 11, 11, 11, 10, 11, 10, 11, 11, 10, 11, 11, 11, 11, 11}

The main disadvantage of stochastic rounding is that the results are not repeatable.

Forms have an extensive set of options to provide fine-grained control over the output. Here I use AccountingForm to display a column of numbers. DigitBlock specifies the grouping factor and NumberPadding allows control of the characters used to pad out the display on the left (shown here as spaces) and right (shown as zeros).

Discussion

Contrast this to AccountingForm without the options.

  In[129]:=  AccountingForm [Column [{100000.00, 1000000.00, 10000000.00}]]
Out[129]//AccountingForm=
             100000.
             1000000.
             10000000.

PaddedForm is convenient when all you want to do is pad out a number with specific characters on the left and right. This is often a useful operation prior to conversion to a string to generate fixed-length identifiers.

Discussion

EngineeringForm forces exponents in multiples of three, provided an exponent of at least three is required.

  In[132]:=  {10.0, 100.0, 1000.0, 10 000.0, 100 000.0, 1 000 000.0} // EngineeringForm
Out[132]//EngineeringForm=
             {10., 100., 1.×103, 10.×103, 100.×103, 1.×106}

ScientificForm always shows numbers with one digit before the decimal and adjusts the exponent accordingly.

  In[133]:=  {10.0, 100.0, 1000.0, 10 000.0, 100 000.0, 1 000 000.0} // ScientificForm
Out[133]//ScientificForm=
             {1.×101, 1.×102, 1.×103, 1.×104, 1.×105, 1.×106}

You can use the option NumberFormat to get precise control of the display. NumberFormat specifies a function (see Chapter 2 for details) that accepts up to three arguments for the mantissa, base, and exponent. Here is an example that displays numbers like a calculator might. Here, the function uses Row to format the mantissa and exponent (it ignores the base).

Discussion