Have you ever noticed that whenever a large prime is found in mathematics, the following sentence always pops up somewhere: “.. and it’s last 10 digits are: xxxxxxxxxx “. How are these terminating digits found? Do they just calculate the whole number and then look at the tail? Does it take a supercomputer to actually find these digits?

Turns out: no, not really. Even you can do it.

Where do the last digits come from?

Let’s first consider where these terminating digits actually come from. Say that we want to find the last 2 digits of 7531 * 8642. Don’t pull out a calculator, instead, think about where these last 2 digits would actually come from. Would they be affected by all digits in our numbers? Of course not, if we visualize this multiplication using the area model, only the red areas affect the last 2 digits of our number:

Area model

As you can see, the last 2 digits of all other products are ’00’. Thus, the product of any two numbers whose last 2 digits are ’42’ and ’31’ will have the same last 2 digits, everytime. Specifically, 42 * 31 = 1302, so the last two digits of 7531 * 8642 are ’02’.

Mersenne Primes

Time to go a little bigger. Mersenne Primes are primes given by \(2^p – 1\) where \(p\) itself is prime as well. Pretty much all of the largest primes we know are Mersenne Primes. In fact, the largest prime we know is a Mersenne Prime: \(2^{57885161} – 1\). Want to know what its last 10 digits are?

The problem is somewhat similar to the one we faced previously, except instead of two numbers, we have to multiply 57885161 numbers together. Luckily, since all of these numbers are the same, there is a quick way to do this! Here’s a Python implementation of the algorithm:

def mypow(x, n):
	if(n == 0):
		return 1
	if(n == 1):
		return x
	if(n == 2):
		return x * x
	if(n % 2 == 0):
		return mypow(x * x, int(n / 2))
	return x * mypow(x * x, int((n - 1) / 2))

However, if we now ask it to calculate mypow(2, 57885161) % 10000000000, nothing happens. That’s because in the end, this query will only return the last 10 digits, but it calculates all of them, which is exactly what we were trying to avoid.

The key now is to split our gigantic power into the string of multiplications it really is, and after a multiplication, drop all digits that aren’t in the last 10. As you saw above, to do this we use % 10000000000 or % 10**10 (\(10^{10}\)). Here’s what the code looks like:

def mypow(x, n, digits):
	if(n == 0):
		return 1
	if(n == 1):
		return x % (10 ** digits)
	if(n == 2):
		return x * x % (10 ** digits)
	if(n % 2 == 0):
		return mypow(x * x, int(n / 2), digits) % (10 ** digits)
	return x * mypow(x * x, int((n - 1) / 2), digits) % (10 ** digits)

As you can see, I went for an extra argument specifying how many digits we are looking for. First, test it a little and see if it works. For example, 2^10 = 1024, so mypow(2, 10, 2) should return 24. After that, try mypow(2, 57885161, 10) - 1 (don’t forget the -1!). It should return 1724285951, the last 10 digits of our prime!

Bigger!

If you have a decent computer, the answer should have appeared in a few seconds. We should be able to calculate even larger Mersenne Primes with this method. Prime(3443958) is the largest known prime for which \(2^p-1\) is prime, but, hypothetically, let’s just assume that prime(10000000) satisfies this condition as well! The number would be astronomically large, WolframAlpha estimates it has approximately \(3.43\cdot{10}^{54012208}\) digits. But in a few seconds, our function tells us the last ten digits of this number are 6819899391!

Does this always work?

You might be wondering whether this works for any sort of large number. Unfortunately, it doesn’t always work out as well as for Mersenne Primes. Generally, any number that can be calculated by multiplications and additions only, is fine. Take, for example, the non-Mersenne but still enormously large prime 19249 * 2^13018586 + 1. The following query gives its last 10 digits: (19249 * mypow(2, 13018586, 10) + 1) % (10 ** 10), returning 8618996737.

Anyway, that’s it for now. Happy digit-hunting.

Meet Stirling’s Formula:

\[ \lim_{n \to \infty} \frac{n!}{\sqrt{2\pi{n}}\cdot n^{n} \cdot e^{-n}} = 1\]

Essentially, this formula says that for very large n, the formula \( \sqrt{2\pi{n}}\cdot{n}^{n}\cdot{e}^{-n} \) approximates \( n! \), since the ratio between these two terms approaches one. This is very useful for calculating large \(n!\) since it is a lot faster.

Faster squaring

But why, you might ask? Calculating \(n!\) takes \(n-1\) multiplications, since you are multiplying n numbers together, while the term \(\sqrt{2\pi{n}}\cdot(\frac{n}{e})^{n}\) takes at least as much, since a) it involves an n-th power, and b) it has a squareroot, and other scary things like \(\pi\) and \(e\).

Turns out, there is a really neat method called exponentiation by squaring that allows us to calculate powers (in our case \((\frac{n}{e})^{n}\) ) a lot faster. Here is how it works:

Let’s say we want to calculate the following power: \(3^8\). If we just multiply eight 3’s together, we’ll need 7 multiplications. However, remember that we can split exponents, specifically:

\[a^{nm} = (a^n)^m\]

We can do the same thing with our power by dividing out the two’s:

\[3^8 = (3^2)^4 = ((3^2)^2)^2\]

This is pretty neat – after all, squaring a number is just one multiplication, so now we only need 3 multiplications rather than 7! Let’s try that again, now with something larger:

\[2^{28} = (2^2)^{14} = ((2^2)^2)^7 = 16^7\]

Now we’ve run into a problem: uneven exponents. We can’t split the 7 into 2 and 3.5 since then we’d have to take roots, and our goal was to find a faster  way to calculate power of n. In fact, since 7 is prime there is no good way to split it into two other integer exponents using the above method. However, not everything is lost! We’ll use this other rule for exponents:

\[a^{n+m} = a^n\cdot{a}^m\]

How about we rewrite 7 to 6 + 1?

\[16^7 = 16\cdot{16}^6 = 16\cdot(16^2)^3 = 16\cdot(256)^3 = 16\cdot{256}\cdot(256)^2 \]

In its final form our equation only needs 3 more multiplications, and it took us another 3 multiplications to get the 16’s and 256’s. That makes 6 multiplications in total, as opposed to the original 27 we’d need for \(2^{28}\). This is pretty cool, now let’s try coming up with a formal algorithm for this method.

The algorithm

Let’s assume we are given a value x and a positive integer exponent n. We’re not going to worry about non-positive exponents since those would be somewhat trivial to handle (Cue “this proof is left as an exercise for the reader.“), and fractional exponents are a whole different beast. In other words, n is positive and whole, or, if you’re a mathematician, \(n\in\mathbb{N}\).

The best way to do this would be by using recursion, since our algorithm continuously splits powers until the exponents are small enough (<= 2). Let’s start by defining a function that does this:

function power(x, n)
	if(n == 1) 
		return x;
	elseif(n == 2) 
		return x * x;
	else ?

Those are the special cases, now we get to the recursion. First, we’ll have to check whether the exponent n is even or uneven. If it is even, we can divide it by two and raise x*x to this power. If it is uneven, we subtract 1 and try again.

function power(x, n)
	if(n == 1) 
		return x;
	elseif(n == 2) 
		return x * x;
	elseif((n % 2) == 0)
		return power(x * x, n / 2);
	else
		return x * power(x, n - 1);

(Where % is the modulo operator, i.e. if c is the result of a % b, then c is rest after division a/b. So (n % 2) == 0 if and only if n is even.)

This actually works fairly well, but we can make a slight optimization. We already noticed how, after dividing our exponent by two, there’s no good way to predict whether the new exponent will be even or uneven. However, when we subtract one we do know. If n is uneven, then n-1 must be even! We can skip a step in our recursion by changing the last line to this:

return x * power(x * x, (n - 1) / 2);

The new exponent (n – 1) / 2 might very well be even as well, but this time there is no simple way to predict that. Anyway, our full (pseudo)code now looks like this:

function power(x, n)
	if(n == 1) 
		return x;
	elseif(n == 2) 
		return x * x;
	elseif((n % 2) == 0)
		return power(x * x, n / 2);
	else
		return x * power(x * x, (n - 1) / 2);

So there you have it! Exponentation by squaring, a simple method to greatly reduce the number of multiplications needed for calculating \(x^n\). If you’re a programmer and you think you should maybe implement this for your projects: don’t worry, most of the time compilers use this (or a similar) method as well. It’s just a nice example of how some creative math can get you a long way when it comes to optimizing code.

Back to Stirling

There are other reasons that Stirling’s Formula (or Stirling’s Approximation) is a useful method to approximate n!, but if you hadn’t noticed yet, it was just an excuse to talk about exponentiation by squaring. I would encourage you to research it on your own, though.