The arithmetic sequence, 1487, 4817, 8147, in which each of the terms increases by 3330, is unusual in two ways: (i) each of the three terms are prime, and, (ii) each of the 4-digit numbers are permutations of one another.

There are no arithmetic sequences made up of three 1-, 2-, or 3-digit primes, exhibiting this property, but there is one other 4-digit increasing sequence.

What 12-digit number do you form by concatenating the three terms in this sequence?

Sage Solution

There are really only a handful of ways to go about this, and the timing isn’t going to be too substantial. So, I’m simply going to write an in-line Sage code block and execute it. The idea is to get a list of primes between 1,000 and 10,000, sort their digit decompositions, and then (after seeing which primes have the same digit decompositions), see which prime subsets have an arithmetic sequence underlying the distances between them. No too horribly interesting, but the code below does the job.

import time
start = time.time()
# make a list of primes
P = prime_range(1000,10000)
# for each prime, make a list of digits, and sort that list
L = [p.digits() for p in P]
for i in range(len(L)):
# make a list of digits that appear for at least 3 primes
M = []
for i in range(len(L)-3):
    if L[i]==L[i+1] and L[i+1]==L[i+2] and L[i+2]==L[i+3]:
        if L[i] not in M: M.append(L[i])
# N is a list of lists of at least 3 primes with the same 4 digits
N = []
for m in M:
for p in P:
    s = p.digits()
    if s in M:
        r = M.index(s)
# Use a Sage combinations iterator to examine sets of 3 primes
for n in N:
    X = combinations_iterator(n, 3)
    for x in X:
        if x[1]-x[0] == x[2]-x[1]: print x
elapsed = time.time() - start
print "time elapsed: %s seconds" % elapsed

When executed, we get the following result.

[1487, 4817, 8147]
[2969, 6299, 9629]
time elapsed: 0.0485339164734 seconds


Many interesting computational problems, such as those on Project Euler require that one find the sum of proper divisors of a given integer. I had a fairly crude brute-force method for doing this, and was subsequently emailed a comment by Bjarki Ágúst Guðmundsson who runs the site He pointed me in the direction of this page and provided some sample code illustrating how such an approach runs asymptotically faster than the approach I had been taking. Awesome! I’m going to expand on that a bit here, providing some mathematical proofs behind the claims and providing code for those who may want to take advantage of this.

The Mathematics Behind It All

Let the function $\sigma(n)$ be the sum of divisors for a positive integer $n$. For example,

It should seem obvious that for any prime $p$ we have $\sigma(p)=1+p$. What about powers of primes? Let $\alpha\in\mathbb{Z}_+$, and then

We’d like to write that in a closed form, i.e., without the “$+\cdots+$”. We use a standard series trick to do that.

\[\begin{align} \sigma(p^\alpha) &= 1+p+p^2+\cdots+p^\alpha\cr p\sigma(p^\alpha) &= p+p^2+p^3+\cdots+p^{\alpha+1}\cr p\sigma(p^\alpha)-\sigma(p^\alpha) &= (p+p^2+\cdots+p^{\alpha+1})-(1+p+\cdots+p^\alpha)\cr p\sigma(p^\alpha)-\sigma(p^\alpha) &= p^{\alpha+1}-1\cr (p-1)\sigma(p^\alpha) &= p^{\alpha+1}-1\cr \sigma(p^\alpha) &=\frac{p^{\alpha+1}-1}{p-1}.\end{align}\]

That solves the problem of finding the sum of divisors for powers of primes. It would be nice if we could show that $\sigma$ is multiplicative on powers of primes, i.e., that $\sigma(p_1^{\alpha_1}p_2^{\alpha_2})=\sigma(p_1^{\alpha_1})\sigma(p_2^{\alpha_2})$. We’ll prove that this is the case, and solve the problem in general along the way.

Proposition: The function $\sigma$ is multiplicative on powers of primes.

Proof: Let $n$ be a positive integer written (uniquely, by the fundamental theorem of arithmetic) as
\[n=\prod_{i=1}^m p_i^{\alpha_i}\]
for $m$ distinct primes $p_i$ with $\alpha_i\in\mathbb{Z}_+$. Any divisor $k$ of $n$ then has the form
\[k=\prod_{i=1}^m p_i^{\beta_i}\]
where each $\beta_i$ satisfies $0\le\beta_i\le\alpha_i$. Then
\[\sigma(n)=\sigma\left(\prod_{i=1}^m p_i^{\alpha_i}\right)\]
is the sum of all divisors $k$ of $n$ and can be written by summing over all possible combinations of the exponents $\beta_i$. There are $\prod_{i=1}^m \alpha_i$ combinations, and we can form their sum and simplify it as follows.
\[\begin{align}\sigma(n) &= \sum_{1\le i\le m,\ 0\le\beta_i\le\alpha_i}p_1^{\beta_i}p_2^{\beta_2}\cdots p_m^{\beta_m}\cr &= \sum_{\beta_1=0}^{\alpha_1}p_1^{\beta_1}\left(\sum_{2\le i\le m,\ 0\le\beta_i\le\alpha_i}p_2^{\beta_2}p_3^{\beta_3}\cdots p_m^{\beta_m}\right)\cr &= \sum_{\beta_1=0}^{\alpha_1}p_1^{\beta_1}\sum_{\beta_2=0}^{\alpha_2}p_2^{\beta_2}\left(\sum_{3\le i\le m,\ 0\le\beta_i\le\alpha_i}p_3^{\beta_3}p_4^{\beta_4}\cdots p_m^{\beta_m}\right) \cr &= \vdots\cr &=\sum_{\beta_1=0}^{\alpha_1}p_1^{\beta_1}\sum_{\beta_2=0}^{\alpha_2}p_2^{\beta_2}\sum_{\beta_3=0}^{\alpha_3}p_3^{\beta_3}\ \cdots\ \sum_{\beta_m=0}^{\alpha_m}p_m^{\beta_m}\cr &=\sigma(p_1^{\alpha_1})\sigma(p_2^{\alpha_2})\cdots\sigma(p_m^{\alpha_m}).\end{align}\]
This completes the proof. Q.E.D.

Thus, we now have a formula for the sum of divisors of an arbitrary positive integer $n$ using the factorization of $n$. Namely,
\[\sigma(n)=\sigma\left(\prod_{i=1}^m p_i^{\alpha_i}\right)=\prod_{i=1}^m\left(\frac{p_i^{\alpha_i+1}-1}{p_i-1}\right).\]

This is something I use quite a bit for various problems and programming exercises, so I figured I could post it here. It’s a basic post that isn’t advanced at all, but that doesn’t mean that the implementation given below won’t save work for others. The idea is to create a list of primes in C by malloc’ing a sieve, then malloc’ing a list of specific length based on that sieve. The resulting list contains all the primes below a given limit (defined in the code). The first member of the list is an integer representing the length of the list.

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#define                 bool    _Bool
static unsigned long    prime_limit = 1000000;
unsigned long sqrtld(unsigned long N) {
    int                 b = 1;
    unsigned long       res,s;
    while(1<<b<N) b+= 1;
    res = 1<<(b/2 + 1);
    for(;;) {
        s = (N/res + res)/2;
        if(s>=res) return res;
        res = s;
unsigned long * make_primes(unsigned long limit) {
    unsigned long      *primes;
    unsigned long       i,j;
    unsigned long       s = sqrtld(prime_limit);
    unsigned long       n = 0;
    bool *sieve = malloc((prime_limit + 1) * sizeof(bool));
    sieve[0] = 0;
    sieve[1] = 0;
    for(i=2; i<=prime_limit; i++) sieve[i] = 1;
    j = 4;
    while(j<=prime_limit) {
        sieve[j] = 0;
        j += 2;
    for(i=3; i<=s; i+=2) {
        if(sieve[i] == 1) {
            j = i * 3;
            while(j<=prime_limit) {
                sieve[j] = 0;
                j += 2 * i;
    for(i=2;i<=prime_limit;i++) if(sieve[i]==1) n += 1;
    primes = malloc((n + 1) * sizeof(unsigned long));
    primes[0] = n;
    j = 1;
    for(i=2;i<=prime_limit;i++) if(sieve[i]==1) {
        primes[j] = i;
    return primes;
int main(void) {
    unsigned long * primes = make_primes(prime_limit);
    printf("There are %ld primes <= %ld\n",primes[0],prime_limit);
    return 0;

Say one wanted to form a list of all primes below 1,000,000. That’s what the above program does by default, since “prime_limit = 1000000.” If one compiles this and executes, you would get something like what follows. The timing is relatively respectable.

$ gcc -O3 -o prime-sieve prime-sieve.c 
$ time ./prime-sieve 
There are 78498 primes <= 1000000
real	0m0.008s
user	0m0.004s
sys	0m0.000s

The code is linked here: prime-sieve.c