Chapter 7: Number theory

Return to all notes

This chapter investigates some aspects of Number Theory, which is the study of the integers. Prime numbers are one of the most important aspects of Number Theory. We will use this topic as an example of scientific computing. We will look at some of the computational study of number theory throughout this chapter, however, in short, this chapter shows a lot of interesting examples and how they are coded. Later we will spend some time investigating the efficiency of the code here.

Prime Numbers

Recall that a prime number is a positive integer whose only factors are 1 and itself. Most languages have built-in versions of testing for primes, and Julia has one built in, however, we will start by writing our own versions.

Let’s first start with a method that returns all of the factors of a given integer:

function all_factors(n::Integer)
  factors = [1]
  for i=2:n-1
    if mod(n,i)==0
      push!(factors,i)
    end
  end
  push!(factors,n)
end

This function starts with an array (list) consisting of only the number 1. It then checks all numbers between 1 and n to see if each is a factor. If it finds a factor, add it to the list with the push function.

Exercise

The template for the function should be

function is_prime(n::Integer)

end

Finding Primes

We’ll now try to find a prime number. Many languages have a function called next_prime that takes a number as an input and returns a prime number greater than that. Here’s a code fragment that may work:

function next_prime(n::Integer)
   i=n+1
   # check if i is prime and if not increment by 1
   # once you find a prime, return it
 end

Exercise

Perfect Numbers

A perfect number is a number, n which has the sum of the factors (except itself) equal to n. For example, 28 is a perfect number because the factors are 1,2,4,7,14,28. The sum of the factors (except itself) is 1+2+4+7+14=28 A nice function would be to check if a number is perfect. Here’s a shell:

function is_perfect(n::Integer)
   # find the factors of n and check the sum
end

Exercise

Happy Numbers

Here’s another fun example. A number \(n\) is called happy by the following
1. Take the digits of \(n\) and square each one.
2. Sum the squares.
3. if the sum is 1, then the number is happy, if not repeat these steps.

Note: it’s also helpful that if this process results in the number 4, then you can never result in a sequence that reaches 1. You can call these number unhappy.

Here’s some examples:

and since we have returned to 4, this will continue cycling, so we stop and say 4 is unhappy. (It has been proven that any unhappy number will eventually hit this cycle.)

To write this in julia, first, it is helpful to have a procedure that splits any positive integer into its digits. Here’s such a procedure:

function split_digits(n::Integer)
  reverse(map(d->parse(Int64,d),split(string(n),"")))
end

Just to understand how it works, do split_digits(1234) which returns the array [4,3,2,1], a array of the ones, tens, hundreds digit, etc.

Next, let’s write a function called is_happy which determines if a number n is happy. The function should be recursive and return true if it is 1, return false if it is 4 and do the steps above otherwise. The template for the function should be

function is_happy(n::Integer)

end

Finding Happy numbers

One way to find happy numbers then is check a range of integers using the function you wrote above. For example, you can type

[is_happy(n) for n=1:100]

which will return an array of true/false if the number is happy or not.

This isn’t so helpful though if you look halfway down the list, what number is happy? A better way to do this is to use julia’s built-in function filter. For example,

filter(is_happy, 1:100)

returns all integers between 1 and 100 that are happy.

Primes can be big!!

If you do a quick google search on the largest prime number known, you’ll see that primes are big. As of early 2016, the largest prime had 22,338,618 digits. Let’s not go that big.

Look at the number \(2^{89}-1\), which is prime. If we type this into julia, we get -1, does this mean that \(2^{89}=0\)? A bit of a hint is to try to raise 2 the 60th power and higher. What happens? Remember that to store a number in a computer, it is typically done using binary. You can use the function bits to see how the number is stored. Try bits(2^60) and higher powers. What’s going on?

Does this mean we’re stuck trying to see if \(2^{89}-1\) is prime. Actually, there is a built-in data type called a BigInt that can handle any sized integer.

Here’s how we can get big integers. Type big(3) to get a big integer with value 3. Type typeof(ans) to show that. If we try big(2^89-1), you will get a big integer with value -1. Why? How can we use big to get that value.

Notice that I haven’t recommended that you use is_prime to test \(2^{89}-1\) or even \(2^{60}-1\). Try timing is_prime(2^{2}8-1) and is_prime(2^{2}9-1). Predict how long \(2^{89}-1\) will take.

We used a very bad method to determine primality and will reexamine it later.

Speeding up is_prime

If you time the is_prime function as we wrote above, you’ll notice how slow it is. The reason why is it relies on the all_factors function which checks every every number between 2 and n if it is a factor. This is very ineffeicient, so we will try some ways to speed things up.

Placing Functions in a Separate file

Note: this isn’t really a topic in Number Theory, but a way to bundle a number of functions together.

Often, there are a number of related functions that you may need. In this case, consider all of the related number theory functions in this chapter. We will put all of these in the same file and load it when needed. Here are the steps:

  1. Open a new text file in JuliaBox or in Jupyter on your local machine. To do this, under NEW (upper right corner of the main page), select Text File. Check under the main menu that your

  2. Change the Language to Julia under the Language pull down menu.

  3. Change the name of the file to filename.jl. where filename is appropriately named.

  4. You can type all of the functions in this file like:

    function all_factors(n::Integer)
    factors = [1]
    for i=2:n
        if mod(n,i)==0
            push!(factors,i)
        end
    end
    factors
end

function is_prime(n::Integer)
    # fill in the details
end

function next_prime(n::Integer)
  # fill in the details
end

function is_perfect(n::Integer)
  # fill in the details
end

  5. To load in these files, type include("filename.jl") in the ipynb file.

  6. Note: if you change anything in the file, you will need to rerun the include command above.