The nth pentagonal number P_n is the number of dots in diagrams like those below with n concentric pentagons.

We have the formula

P_n = (3n² − n)/2

where n is a positive integer. If n is an integer but not positive, the equation above defines a generalized pentagonal number.

If you’re given an n, you can easily compute P_n. But suppose you’re given a large number x. How would you determine if it is a pentagonal number? And if it is a pentagonal number, how would you find n such that x = P_n?

If

x = P_n = (3n² − n)/2

then we can solve a quadratic equation for n:

n = (1 ± √(24x + 1))/6.

If 24x + 1 is not a perfect square, n is not an integer and x is not a pentagonal number, ordinary or generalized. For example,

√(24 × 20260615 + 1)) = 22051.185…

and so 20260615 is not a pentagonal number nor a generalized pentagonal number.

Now suppose

x = 170141183460469231731687303715884105727.

Is this a pentagonal number? You can’t just compute √(24x + 1) in floating point arithmetic because the result is a 20-digit number, and floating point number have 15 digits of precision, so you can’t tell whether the result is an integer.

However, you can compute

⌊√(24x + 1)⌋

with only integer arithmetic using the sqrt_floor function from this post.

def sqrt_floor(n):
    a = n
    b = (n + 1) // 2
    while b < a:
        a = b
        b = (a*a + n) // (2*a)
    return a

The following prints a positive number,

x = 2**127 - 1
y = 24*x + 1
r = sqrt_floor(y)
print(y - r**2)

which tells us y is not a perfect square.

Now suppose y is a perfect square. Then

(1 ± √(24x + 1))/6

is rational, does it have to be an integer? In fact it one, and only one, of the roots will be an integer. If

24x + 1 = r²

then r is congruent to ±1 mod 6 because the left side is congruent to 1 mod 6. If r = 1 mod 6 then the smaller root is an integer, and if r = 5 mod 6 then the larger root is an integer.

So if 24x + 1 = r², then x is a pentagonal number if r = 5 mod 6 and a generalized pentagonal number otherwise.

The function pentagonal_index takes a number x and return n if x = P_n and None if no such n exists.

def pentagonal_index(x):
    y = 24*x + 1
    r = sqrt_floor(y)
    if r*r != y:
        return None
    if r % 6 == 5:
        return (1 + r) // 6
    else:
        return (1 - r) // 6

We can test this with the following code.

P = lambda n: (3*n**2 - n) // 2
for n in [2, 3, -4, -5, 10**200]:
    assert(pentagonal_index(P(n)) == n)

Note that P(10**200) is too big to fit into a float, but the code works fine. This is because we use integer division (//) everywhere. If we had said

P = lambda n: (3*n**2 - n) / 2

the test above would pass for the small values of n but output

OverflowError: integer division result too large for a float

when it came to n = 10²⁰⁰.

Related posts

• Partitions and pentagons
• Certified Fibonacci numbers
• Truncated triangular numbers

There’s a lot of mathematics just off the beaten path. You can spend a career in math and yet not know all there is to know about even the most basic areas of math. For example, this post will demonstrate something you may not have seen about decimal forms of fractions.

Let p > 5 be a prime number and 0 < k < p. Then the digits in k/p might be the same for all k, varying only by cyclic permutations. This is the case, for example, when p = 7 or p = 17. More on these kinds of fractions here.

The digits in k/p repeat for every k, but different values of k might have sequences of digits that vary by more than cyclic permutations. For example, let’s look at the values of k/13.

>>> for i in range(1, 13):
...   print(i/13)
...
 1 0.0769230769230769
 2 0.1538461538461538
 3 0.2307692307692307
 4 0.3076923076923077
 5 0.3846153846153846
 6 0.4615384615384615
 7 0.5384615384615384
 8 0.6153846153846154
 9 0.6923076923076923
10 0.7692307692307693
11 0.8461538461538461
12 0.9230769230769231

One cycle goes through the digits 076923. You’ll see this when k = 1, 3, 4, 9, 10, or 11. The other cycle goes through 153846 for the rest of the values of k. The cycles 076923 and 153846 are called the distinct repeating sets of 13 in [1].

If we look at fractions with denominator 41, thee are six distinct repeating sets.

02439
04878
07317
09756
12195
14634
26829
36585

You could find these by modifying the Python code above. However, in general you’ll need more than default precision to see the full periods. You might want to shift over to bc, for example.

When you look at all the distinct repeating sets of a prime number, all digits appear almost the same number of times. Some digits may appear one more time than others, but that’s as uneven as you can get. A corollary in [1] states that if p = 10q + r, with 0 < r < 10, then 11 − r digits appear q times, and r − 1 digits appear q + 1 times.

Looking back at the example with p = 13, we have q = 1 and r = 3. The corollary says we should expect 8 digits to appear once and 2 digits to appear twice. And that’s what we see: in the sets 076923 and 153846 we have 3 and 6 repeated twice and the remaining 8 digits appear once.

In the example with p = 41, we have q = 4 and r = 1. So we expect all 10 digits to appear 4 times, which is the case.

Related posts

• Cyclic fractions
• Periods of fractions
• Calendars and continued fractions

[1] James K. Schiller. A Theorem in the Decimal Representation of Rationals. The American Mathematical Monthly
Vol. 66, No. 9 (Nov., 1959), pp. 797-798

I recently found out about Andrica’s conjecture: the square roots of consecutive primes are less than 1 apart.

In symbols, Andrica’s conjecture says that if p_n and p_n+1 are consecutive prime numbers, then

√p_n+1 − √p_n < 1.

This has been empirically verified for primes up to 2 × 10¹⁹.

If the conjecture is true, it puts an upper bound on how long you’d have to search to find the next prime:

p_n+1 < 1 + 2√p_n  + p_n,

which would be an improvement on the Bertrand-Chebyshev theorem that says

p_n+1 < 2p_n.

Pentagonal numbers are truncated triangular numbers. You can take the diagram that illustrates the nth pentagonal number and warp it into the base of the image that illustrates the (2n − 1)st triangular number. If you added a diagram for the (n − 1)st triangular number to the bottom of the image on the right, you’d have a diagram for the (2n − 1)st triangular number.

In short,

P_n = T_{2n − 1} − T_{n − 1}.

This is trivial to prove algebraically, though the visual proof above is more interesting.

The proof follows immediately from the definition of pentagonal numbers

P_n = (3n² − n)/2

and triangular numbers

T_n = (n² + n)/2.

Related posts

• Partitions and pentagons
• Tetrahedral numbers
• A pentagon, hexagon, and decagon walk into a bar …
  

If a number has a finite but nonzero fractional part, so do all its powers. I recently ran across a proof in [1] that is shorter than I expected.

Theorem: Suppose r is a real number that is not an integer, and the decimal part of r terminates. Then r^k is not an integer for any positive integer k.

Proof: The number r can be written as a reduced fraction a / 10^m for some positive m. If s = r^k were an integer, then

10^mk s = a^k.

Now the left side of this equation is divisible by 10 but the right side is not, and so s = r^k must not be an integer. QED.

The only thing special about base 10 is that we most easily think in terms of base 10, but you could replace 10 with any other base.

What about repeating decimals, like 1/7 = 0.142857142857…? They’re only repeating decimals in our chosen base. Pick the right base, i.e. 7 in this case, and they terminate. So the theorem above extends to repeating decimals.

[1] Eli Leher. √2 is Not 1.41421356237 or Anything of the Sort. The American Mathematical Monthly, Vol. 125, No. 4 (APRIL 2018), page 346.

The last few posts have looked at expressing an odd prime p as a sum of two squares. This is possible if and only if p is of the form 4k + 1. I illustrated an algorithm for finding the squares with p = 2²⁵⁵ − 19, a prime that is used in cryptography. It is being used in bringing this page to you if the TLS connection between my server and your browser is uses Curve25519 or Ed25519.

World records

I thought about illustrating the algorithm with a larger prime too, such as a world record. But then I realized all the latest record primes have been of the form 4k + 3 and so cannot be written as a sum of squares. Why is p mod 4 equal to 3 for all the records? Are more primes congruent to 3 than to 1 mod 4? The answer to that question is subtle; more on that shortly.

More record primes are congruent to 3 mod 4 because Mersenne primes are easier to find, and that’s because there’s an algorithm, the Lucas-Lehmer test, that can test whether a Mersenne number is prime more efficiently than testing general numbers. Lucas developed his test in 1878 and Lehmer refined it in 1930.

Since the time Lucas first developed his test, the largest known prime has always been a Mersenne prime, with exceptions in 1951 and in 1989.

Chebyshev bias

So, are more primes congruent to 3 mod 4 than are congruent to 1 mod 4?

Define the function f(n) to be the ratio of the number of primes in each residue class.

f(n) = (# primes p < n with p = 3 mod 4) / (# primes p < n with p = 1 mod 4)

As n goes to infinity, the function f(n) converges to 1. So in that sense the number of primes in each category are equal.

If we look at the difference rather than the ratio we get a more subtle story. Define the lead function to be how much the count of primes equal to 3 mod 4 leads the number of primes equal to 1 mod 4.

g(n) = (# primes p < n with p = 3 mod 4) − (# primes p < n with p = 1 mod 4)

For any n, f(n) > 1 if and only if g(n) > 0. However, as n goes to infinity the function g(n) does not converge. It oscillates between positive and negative infinitely often. But g(n) is positive for long stretches. This phenomena is known as Chebyshev bias.

Visualizing the lead function

We can calculate the lead function at primes with the following code.

from numpy import zeros
from sympy import primepi, primerange

N = 1_000_000
leads = zeros(primepi(N) + 1)
for index, prime in enumerate(primerange(2, N), start=1):
    leads[index] = leads[index - 1] + prime % 4 - 2

Here is a list of the primes at which the lead function is zero, i.e. when it changes sign.

[   0,     1,     3,     7,    13,    89,  2943,  2945,  2947,
 2949,  2951,  2953, 50371, 50375, 50377, 50379, 50381, 50393,
50413, 50423, 50425, 50427, 50429, 50431, 50433, 50435, 50437,
50439, 50445, 50449, 50451, 50503, 50507, 50515, 50517, 50821,
50843, 50853, 50855, 50857, 50859, 50861, 50865, 50893, 50899,
50901, 50903, 50905, 50907, 50909, 50911, 50913, 50915, 50917,
50919, 50921, 50927, 50929, 51119, 51121, 51123, 51127, 51151,
51155, 51157, 51159, 51161, 51163, 51177, 51185, 51187, 51189,
51195, 51227, 51261, 51263, 51285, 51287, 51289, 51291, 51293,
51297, 51299, 51319, 51321, 51389, 51391, 51395, 51397, 51505,
51535, 51537, 51543, 51547, 51551, 51553, 51557, 51559, 51567,
51573, 51575, 51577, 51595, 51599, 51607, 51609, 51611, 51615,
51617, 51619, 51621, 51623, 51627]

This is OEIS sequence A038691.

Because the lead function changes more often in some regions than others, it’s best to plot the function over multiple ranges.

The lead function is more often positive than negative. And yet it is zero infinitely often. So while the count of primes with remainder 3 mod 4 is usually ahead, the counts equal out infinitely often.

If p is an odd prime, there is a theorem that says

x² = −1 mod p

has a solution if and only if p = 1 mod 4. When a solution x exists, how do you find it?

The previous two posts have discussed Stan Wagon’s algorithm for expressing an odd prime p as a sum of two squares. This is possible if and only if p = 1 mod 4, the same condition on p for −1 to have a square root.

In the previous post we discussed how to find c such that c does not have a square root mod p. This is most of the work for finding a square root of −1. Once you have c, set

x = c^k

where p = 4k + 1.

For example, let’s find a square root of −1 mod p where p = 2²⁵⁵ − 19. We found in the previous post that c = 2 is a non-residue for this value of p.

>>> p = 2**255 - 19
>>> k = p // 4
>>> x = pow(2, k, p)

Let’s view x and verify that it is a solution.

>>> x
19681161376707505956807079304988542015446066515923890162744021073123829784752
>>> (x**2 + 1) % p
0

Sometimes you’ll see a square root of −1 mod p written as i. This makes sense in context, but it’s a little jarring at first since here i is an integer, not a complex number.

The next post completes this series, giving a full implementation of Wagon’s algorithm.