Tutorial

The tutorial has several sections. You can open a section by double-clicking on the
 arrow at the lower right corner of the section. Or, you can select a section and use the 
 command open all subgroups from the Cell menu.

Introduction

We can use the computer algebra system Mathematica  to investigate and solve problems in 
 number theory. First, we describe some basic features of Mathematica and then
 more specialized aspects that are particularly useful for number theory. After this brief
 tutorial, we give a detailed description of list manipulation, functions, and procedural
 programming.  It is best to work through all the examples yourself.

  Since you are reading  the notebook version of this manual, you already know how
 to start Mathematica.  A command in Mathematica is evaluated by using the
 Enter key  (Macintosh systems) or Shift + Return on Windows and Unix front
 ends. If you wish to run  Mathematica from the command line, then a statement
 can be  evaluated as soon as it is complete by pressing the return key.

Mathematica as a Calculator

The simplest use of  Mathematica is as a calculator.  Unlike standard calculators, there
 is no limit to the size of numbers or the type of functions that Mathematica can handle.
 An expression is evaluated by entering it in the usual mathematical notation. You can evaluate
 the following expression by pressing the shift and return keys simultaneously.

2+ 11 * 8

Try other expressions to get used to the  syntax used in Mathematica.   Here are 
 a few examples.

2^67 -1

(11^23-6^23)/5

17/29

Notice that  the result of 17/29 is again 17/29. This is an important feature of Mathematica.
Mathematica represents mathematical values exactly,  without any approximations. The 
 value of 17/29 is 17/29, so that is the result given by Mathematica. If you wish to see
 a decimal approximation, then you can use the N function of Mathematica.  The
 command
  N[ expression, places] gives a approximation of the expression to the specified number
 of places.

N[17/29, 8]

Exercise: Compute the value of square root of 2 to 50
 decimal places.  The Mathematica function for evaluating the square root is
 Sqrt[  expression].

Mathematical constants like e, pi, and i are represented
 by E, Pi, I in Mathematica.  Mathematica has all built-in variables and functions starting with a 
 capital letter, hence it is a good idea to start all your variables and functions with a lower
 case letter to avoid any confusion with an existing name.

Another feature of Mathematica is worth  paying  attention.  Recall that we use the symbol
 ^ for exponentiation. Let us compute the cube root of -1.

(-1)^(1/3)

N[ (-1)^(1/3), 6]

Observe that the approximation of the cube root of -1 is a complex number. The number
-1 has three cube roots, two of them complex. Mathematica works internally with complex
 numbers, so  it has selected  a primitive complex cube root.  To be able to compute
 a real cube root of a  negative real number, we have to define a function ourselves. We 
 will see how to accomplish this in a later section.

Exercise: Determine sin(pi/12) to 8 decimal places. Compute sin(x)/x for values of
 x close to zero to observe the limit as x approaches 0.  The Mathematica function for
evaluating  sin(x) is Sin[x].

Here are examples of a few more useful functions.  Quotient[a,b] gives the quotient
 when a is divided by b and Mod[a,b] gives the remainder that has the same sign as b. The 
 two functions satisfy the relation

a= b Quotient[a,b] + Mod[a,b].

Quotient[ -34, 7]

Mod[ -34, 7]

Quotient[ 34, -7]

Mod[ 34, -7]

Exercise: Verify that the relation between   the Quotient and  Remainder that is
 specified above is valid.

Another useful function is  Floor[ x], that returns the  largest integer less than x.  Here x has to
 be an expression that evaluates to a real number, otherwise Floor doesn't do anything.
Related functions are Ceiling and Round; Ceiling returns the smallest integer  greater 
 than the expression and Round returns the nearest integer.  Here are a few examples.

Floor[ -4.5]

Ceiling[-4.5]

Round[ -4.6]

Round[ -4.4]

Exercise: What does Round do for  numbers of the form  integer + 1/2?

Floor[1/3]

Floor[Sqrt[2]]

Brackets, Parenthesis, and Braces

Brackets, parenthesis, and braces are intended for  different purposes. Each has
 a specific meaning and their use cannot be interchanged.

Paranthesis:   These are used for grouping expressions. Careful use of parenthesis
 is necessary to make clear the meaning of arithmetical expression.  This is because multiplication
 and division have a higher order of precedence than addition and subtraction.
 For example, 3 + 7* 8 is not (3+7) *8 but  3 + (7*8). The expression 
 21/7-5 is  (21/7) -5 and not  21/(7-5).

21/7-5

21/(7-5)

In Mathematica (as in other programming languages) it is essential to use parenthesis
  to make the meaning of arithemetical expressions transparent.

Exercise: What does  a/b/c represent in Mathematica; (a/b)/c or a/(b/c)? What about
 a/b*c?

Brackets:  Square brackets are used for specifying arguments of functions. For example,
 Floor takes a single argument, Floor[ x]. Square brackets cannot be used for grouping terms
 in arithmetical functions as done in written mathematical work.  Try it in a few examples to 
 see what happens.

2*[3+4]

Braces:   These are used for specifying lists.  A list is an order collection of objects. Lists are
 used for representing vectors, matrices  and other mathematical collections.  Many
 functions will return lists as their result and other functions manipulate lists. List
 manipulation is one of the most powerful features of Mathematica and is discussed in more
 detail in a later section. Here are a few examples of their use.
 We can compute the remainder obtained by  the first few Fibonacci numbers  using the list 
 manipulation features of Mathematica.

numbers= { 1,1,2,3,5,8,13,21,34};

Mod[ numbers, 3]

Notice the semicolon at the end of the definition of the list named numbers above. The
 purpose of the semicolon is suppress the output of the expression.  Try this in a few of
 your own computations.

Referring to previous results

Mathematica  has a convenient mechanism to refer to previously computed
 results in a given session. The symbol % can be used to refer to  the last computed
 result. %n can be used to access the result  in the cell Out[n].

Prime[ 5]

PrimePi[ %]

Here, Prime[n] returns the nth prime and PrimePi[x] returns the number of
 primes less than or equal to x.

 The use of % and %n in expressions should be used very carefully since % refers to
 the last computed output and not the one right above the  current expression. If you
 do another computation between using % and output you want to refer to, then you will
 get incorrect results.

Getting Help

Mathematica has extensive online help to access information about functions
 and  commands. The ? operator can be used for getting help. If a function 
 is prefixed with ?, then Mathematica prints usage information.

?Quotient

More detailed information can be obtained by using ??.

??Mod

If you don't know the name of a command, then you can
 search by using a * in the command name. For example,
 to find all commands beginning with Div, use ?Div*

?Div*

Exercise: Find out what DivisorSigma does and use to verify that 2305843008139952128
 is a perfect number.  You can also use Divisors to find the divisors of the number and
 FactorInteger to find the prime factorization.

Functions

The power of Mathematica  lies in the enormous number of built-in functions and
 the numerous  features  that can be combined to define  complex new functions that
 significantly extend the computational capabilities of Mathematica. This allows the 
 user to quickly  expand the abilities  of Mathematica in powerful ways. 

 A function is defined by specifying the name and its argument list, followed by  actions
 to be performed by the function. The syntax  will be clear from a few examples.

square[x_]:= x^2;

square[ 5]

square[ Sqrt[2]]

In the left side of the definition, we specify the arguments of the function by
 placing an underbar after the variable. The := specifies that it is the definition of
 a function. After this, we give the rule for the evaluation of a function.

Exercise: Define a function to cube of a number.

It is essential to use the underbar to  represent a variable. If you  use an underbar,
 then any number can be used as an argument. Otherwise,  only the argument will
 match the  expression. For example,

squareroot[x]:= x^(1/2)

squareroot[4]

squareroot[x]

In this case, only when the argument of squareroot is x does the function
 give the correct value. No other expression matches it because of the lack of
 the underbar. 

The simplest way to extend the function capabilities is to  able to define functions
 with different expressions for different ranges of the input variables. Recall that we
 did not get the real cube root of -1 in a previous computation.  Since the value of
  the cube root of a negative number is the negative of the cube root of its absolute 
 value, we can define a function to give the real valued cube root as follows.

realcuberoot[x_]:= (x)^(1/3) /; x>=0;
realcuberoot[x_]:= - (-x)^(1/3) /; x<0;

This defines a real valued cube root for real numbers. Notice, the definition for 
 negative numbers. The /; after the expression specifies the conditions under which t
 the function is to be evaluated. Let us try this on a few values.

realcuberoot[-1]

realcuberoot[ -8]

realcuberoot[ 27]

We can  define functions by specifying rules for evaluating the functions. To define
 the Fibonacci numbers, we can convert  the standard  recurrence relation that defines
 them  into a Mathematica function.

fib[1]=1;
fib[2]=1;
fib[ n_]:= fib[ n-1] + fib[n-2]

fib[10]

fib[30]

If you evaluated fib[30] then you would  have noticed that the computation took a few
 minutes (depending on the speed of your machine). This happens because to evaluate
 fib[30], Mathematica first computes fib[29] and fib[28]. To evaluate fib[29] it
 has to compute fib[28] and all the smaller Fibonacci numbers, but these  intermediate
 computations are then thrown away. To compute fib[28], then Mathematica again
 computes all the previous values. This is very inefficient in terms of computational time,
 so there is a way to tell Mathematica to remember all the  intermediate values in
 a computation.  This is done  in the following way. We first  make Mathematica
 forget the old definition by using the Clear command.

Clear[fib]
fib[1]=1;
fib[2]=1;
fib[n_]:= fib[n]= fib[n-1] + fib[n-2]

Here  the definition fib[n_]:= fib[n] tells Mathematica to remember all the
 intermediate values. Now the computation should work much faster.  There are also
 ways to compute large Fibonacci numbers directly, without computing all the
 previous ones.

Exercise:  Define functions to compute n! and |x|.

More detailed description of functions using  the Module command are
 described in the section on procedural programming.

Number Theoretic Functions

Some of the  built-in functions that are  useful for doing number theoretic work
 are the following. 

Mod[k,n] 	: computes the remainder when k is divided by n.

Quotient[k,n]: Gives the integer quotient defined by Floor[k/n].

GCD[n1,n2,...]:  Gives the greatest common divisor of the arguments.

PrimeQ[n]:     Returns True if the number is prime and False otherwise.

Prime[k]:       Gives the kth prime.

FactorInteger[ n]: Returns the prime factorization of n. The function is only suitable
                          for integers less than 10^12. 

Other useful functions are PowerMod and EulerPhi.

Here are a few examples of using these functions.

PrimeQ[ 2^61 -1]

This shows that  2^61 -1 is prime.  Use of PrimeQ is not a proof  that a number is
 prime.   For proofs of primality of a number, consult your number theory textbook.

FactorInteger[ 1234567]

The result of   FactorInteger is a list of the prime factors and the exponent to which
 that prime occurs in the factorization. In the above, 1234567 = 127 * 9721.

 FactorInteger is limited to about 12-15 digits and for numbers with more digits other
 methods are required. Consult  your number theory book for details on other methods.
 Many functions such as Divisors, DivisorSigma and EulerPhi depend on FactorInteger,
 hence they too are restricted in the size of the input.

For small numbers, we can verify if a number n is prime by dividing it
 by all  primes up to sqrt(n).  This is a good illustration of combining Mathematica functions.
 Let us verify that 9721 is Prime. We compute the primes up to sqrt(9721) using the following c
commands.

n=9721;
k=N[ Sqrt[9721]]

numprimes= PrimePi[ k]

listofprimes= Table[ Prime[ i], { i, 1, numprimes}]

Next, we can divide 9721 by all these numbers. Many Mathematica functions
 are Listable,  that is, they can take a list as an argument and apply the function
 to each element. ( The same thing would require the use of  loop structures in
 other languages.) The Table command and  few other list related functions are discussed in
 the following section.

Mod[ n, listofprimes]

We  computed the remainder when n is divided by all the primes up to square root of
 n and we see that all the remainders are non-zero.  This proves that the number  9721
 is prime.

Exercise: Verify that 7927 and  47419 are primes by using this division method.

This method is not suitable for large numbers.  More efficient methods for large primes
 are discussed in your textbook.

Lists

Mathematica  is a powerful list processing language and has a lot of commands 
 to manipulate and create lists.  A list is an ordered collection of objects and Mathematica
  has commands for inserting elements to a list, to manipulate and apply functions to
 the list. 

 The simplest list creation function is the Table commad.  The first argument of Table
 is the function and the second is a range. To create a table of  the first  ten squares, we can
 write

Table[ x^2, { x, 1, 10}]

The range  { x, 1, 10} specifies that x ranges from 1 through 10,  increasing by 1 in
 each step. If we wanted the first ten odd squares, then we can use the range { x, 1, 20, 2}.
 Now  x is incremented by 2 in each step. An ith element of a list L can be accessed by
  L[[i]]. The first and last elements of L are obtained by First[L] and Last[L]. A list is like a
 vector and a list of lists (all identical size) is like a matrix, so matrix operations
 can be applied to these.  Transpose  exchanges the rows and the columns.

Consider the application of these in the following example.

factors= FactorInteger[ 789672]


x= Transpose[factors]

primefactors= First[x]

The result of  these operations is to extract the prime factors of the number. FactorInteger
 returns a list of lists, ( a matrix) whose rows have two elements each, a prime and its
exponent in the factorization.  Transpose interchanges the rows and columns, so the first
 element of the transposed list is now a list of the  prime factors. 

 More list processing functions are discussed in a separate chapter.

Exercise:  Make a list of the first fifty Fibonacci numbers. Divide them by3.
 Which are divisible by 3? Make a conjecture and test it for larger Fibonacci numbers
 by  making a list of only those Fibonacci numbers that you think are multiples of 3. 
 Repeat the exercise to test divisibility by 3.

The Clear command

An important feature of  Mathematica is that it  remembers the values of variables
 and definitions of functions throughout the session. This can be a problem when the
 same names are used again,  because Mathematica will use the old values.  You can
 remove the values and names from Mathematica's memory by using the Clear command.
It is important to do this, as Mathematica will not compute a new value for a function
 that is stored in memory.  Not clearing names is the most frequent source of erroneous
 results in Mathematica programs. Consider the following example.

Suppose  f is a function for which f[2]=4 is already defined.

Let us redefine f and compute f[3];

f[3]=8;
f[x_]:=x^3;
f[3]

Let us  redefine f to a  new function and see what happens.

f[a_]:= a+2;
f[3]

f[2]

f[5]

The redefinition does not effect the values of the function. For this, it is necessary to clear
 all the values if the function. It is best to place the clear statement before the definition
 of a function. This way, if you make changes to the function, then  the function is cleared
 before the changes are entered and all values will be computed again.

Clear[f]
f[a_]:= a+2;
f[3]

Loops and Conditionals

If statement

The If statement is used to perform actions based on the results of tests. The general
 format of the statement is 

 If[ test, statements1, statements2]

 where statesments1 is performed if the  result of the test is True and statements2
 is performed if the result of the  test is False. If nothing is to be done when the test
 is False  then statements2 can be omitted. It is important to separate the statements
 with a comma. The statements can be a series of expression separated by semi-colons.

Clear[ realcuberoot]
realcuberoot[x_]:= If[ x>= 0,
                         x^(1/3),    (* True value *)
                         -(-x)^(1/3) (*False value *)
                      ]

realcuberoot[-1]

The function returns the real cube root of a negative number. Observe that we
 placed a comment about the result in the If statement. Mathematica comments
 are placed between the expressions (*  and *). 

 Another example of an If statement  is the following in which we check if a number
 is a perfect square.

squareQ[ n_]:= If[ Floor[ N[Sqrt[n]]]^2== n,
                        True,  (* Yes, aperfect square *)
                        False (* not a square *)
                  ]

squareQ[196]

squareQ[255]

The function can be written without any If statements by writing
  squareQ[n_]:= Floor[ N[ Sqrt[n]]]^2 == n;

 This is because the result of the test for equality is a True or  False. The following are
 the relational operators that can be used. These can be combined with logical operators,
 And, Or , or Not to create all the tests that we need.

 x = = y       Is x equal to y?

x!= y           Is x not equal to y?

 x> y        Is x greater than y?

x< y          Is x  less than y.

 Similarly, the relations >= and <= are self-explanatory.  The only issue to keep in mind is
 that  two equal signs are necessary for a test of equality. A single = is an assignment of
 one to the other and not a test.  The boolean AND operator is represented by
  && and the logical OR by ||. The complement is represented by the exclamation !.
 For example, to test  if n is a prime that is not of the form 5k+3, we can use

 PrimeQ[n]  &&  !( Mod[ n,5]==3).

Exercise: Write a function that returns True if a number is either a prime or a perfect square, 
 and False otherwise..

Write another function that returns True if a number is a prime  or the square of a
 prime, and False otherwise.

The While statement

The While  statement  is used to evaluate a series of statements as long as
 the specified condition is True. The  structure of the While statement is the following.

 While[ condition,
            statements]

 The statements are repeatedly executed while the condition holds True and the
 loop is exited when the condition is False. 

 Consider the following function to find the first prime larger than an integer. We use 
PrimeQ to test if a number is prime. The function also illustrates the use of the 
 Module function to define  functions.

firstprimeabove[ n_]:= Module[ {k},
                              k=n+1;
    					      While[ !PrimeQ[k], 
                                     k=k+1;
                                   ];
                              Return[k];
                              ]

firstprimeabove[ 167]

In the While  loop, we increment k  as along as it is not prime. When the loop is 
 finished k is prime.  Observe that even if  the input is prime, then the function returns the
 next  prime. The first statement of Module is a list of the variables internal to the Module, and
 are placed within braces. After that, there is a series of statements separated by semi-
colons.

Exercise:   Modify the function to return the second prime above  n,  or more
 generally, the function should take a second argument r, and return the rth prime
 above n.

Do and For loops:

The Do loop  executes  a series of statements  for a  variable in the specified range.

The general format is 

 Do[ statements,
          range]

 As an example, suppose we wish to print all the primes less than 100 that are
 of the form 8k+1, then we write

Do[ If[ PrimeQ[n] && (Mod[n,8]==1),
              Print[ n]
       ],
     {n, 3, 100,2}]

Exercise: The condition in the If statement can be reduced to PrimeQ[n] by
 modifying the range of n. Explain how to do this.

Exercise:   Write a program to print primes of
 the form x^2 +y^2. You can use one Do loop within another.

The For Statement

The For  statement is more general than the Do loop and equivalent to the While
 statement.   As a generalization of the Do it also includes a test and an increment.

 The format is
  For[ init, test, increment, statements]
 where init is the initialization sequence,  and then the statements are executed while
 the test is true. The increment is performed after the statements are executed and the
 test is performed before the statements are executed.

For[ x=1, x< 10, x++, Print[ Prime[x]] ]

Here  we used the increment operator ++ to increase the value of x by one. x++ is
 equivalent to x= x+1;

The Print statement

We have already seen a few examples of the use  of the Print statement. As Mathematica
 prints the result of most calculations, it is usually not necessary to use the Print statement, except
 to format the output  for improved readability.  There are some features of the
 Print statement that should be kept in mind. The first is that Print does not leave any
 space between the values of different arguments. The amount of space for proper
 display has to be specified by you. For example,

Print[ a,b,c]

To get  the appropriate space, specify it in the statement.

Print[ a, " ", b, " ", c]

Additional characters that are useful are "\t" for the tab and "\n" for the end of line. Unlike
 C, Mathematica always places an end of line character at the end of a Print statement.

Do[ Print[ "First prime above ", 10^n, " is ",
              firstprimeabove[ 10^n]],
     {n, 1, 10}]

File IO

Sometimes it is necessary to write data to a file or to read data from a file. An 
 expression can be stored in a file using the  redirection arrows, >>. For example,

{ 2, 3, 5, 7, 11, 13, 17, 19,23}>> primefile

will create  a file called primefile where the list will be stored. This can be  read into
 any expression using the ReadList command.

primes= ReadList[ "primefile", Expression]

Notice that this has  placed an extra set of braces around the input. The extra
 set of braces can be removed using the Flatten command, which removes all the
braces (except the outermost ones) from a list.

primes= Flatten[ primes]

To add to a file  we can use the PutAppend command or its short form, >>>.

{ 29, 31, 37, 41, 43, 47} >>> primefile

primes= Flatten[ ReadList[ "primefile", Expression]]