There’s a problem that a coworker recently brought up. The problem statement can be found here.

This is a re-hash of an existing problem known as the one hundred prisoners problem. It seems impossible at first, and many naive approaches yield failing results. But, given the correct approach, it is solvable with a surprisingly high probability.

A Solution

It’s really just a permutation problem. The basic idea is that the boxes are numbered 1..100 and they all contain a bill with a number 1..100. To construct a working solution, you use the bills as pointers to box numbers in the following way. You start at the box labelled with your number, and go to the box number of the bill it contains. Keep doing that, looking in the box and going to the box of the bill you find inside. Eventually, you get to the box containing the bill that points to where you started – and that’s the bill you’re looking for. So, the question really becomes: How often does a random permutation of degree $n$ contain a cycle having length greater than $n/2$? If you get a longer cycle, then anyone with a number in that cycle is screwed because it takes them longer than $n/2$ steps to find their bill. If there are only cycles having lengths less than $n/2$, then everyone can find their bill in the appropriate number of steps.

So, how many random permutations of degree $n=100$ contain only cycles having length less than 50? For this, we can look at approximations given by Stirling numbers of the first kind. Essentially, the number of random permutations having only cycles of length less than $n/2$ (in our case 50) as asymptotic to $log(2)$, which is around 0.3010. That’s a bit surprising, because it means that the strategy above is expected to win roughly 30% of the time for large numbers of boxes.

Testing with Sage

How well does the strategy actually perform when $n=100$? We can test this out in Sage.

def max_length_orbit_random_perm(n):
    """ Generate a random permutation on 1..n and compute the longest orbit."""
    return max([len(i) for i in Permutations(n).random_element().cycle_tuples()]) 
 
won = 0
tries = 10000
for i in range(tries):
    if max_length_orbit_random_perm(100) <= 50:
        won += 1
print "won %s out of %s" % (won, tries)

And we get the following.

won 3072 out of 10000

So, we win roughly 3 out of 10 times. If we were charged $1 to play the game, but won $100 if everyone found their dollar, then we’d be doing very well to play the game.


Problem

An irrational decimal fraction is created by concatenating the positive integers:

\[0.123456789101112131415161718192021...\]

It can be seen that the 12th digit of the fractional part is 1.

If $d_n$ represents the nth digit of the fractional part, find the value of the following expression.

\[d_1\times d_{10}\times d_{100}\times d_{1000}\times d_{10000}\times d_{100000}\times d_{1000000}\]

Python Solution

This was coded in Python using the Sage Notebook (for speed). The idea is to concatenate the string representations of increasing integers into a single string, turn that string into a list with the appropriate index range, and take the appropriate product from within that list. I use the range 400000 to form the string, guessing (correctly) that this would provide enough digits in the resulting string/list.

import time
 
start = time.time()
 
s = ""
for i in range(400000):
    s += str(i)
d = list(s)
n = prod([ZZ(d[10**i]) for i in range(7)])
 
elapsed = time.time() - start
 
print "result %s found in %s seconds" % (n, elapsed)

When executed, we get the following result.

result 210 found in 0.369875907898 seconds

This is just another example of a problem I decided to attempt as a break at work, trying to form the solution before the cookie on the Project Euler site ran out.

Problem

It can be seen that the number, 125874, and its double, 251748, contain exactly the same digits, but in a different order.

Find the smallest positive integer, x, such that 2x, 3x, 4x, 5x, and 6x, contain the same digits.

Python Solution

This could be written in C or Cython to perform much more quickly, and I’m sure there are other optimizations that could be made to this code, but this was written and executed in a total of about three minutes using the Sage Notebook.

import time
 
start = time.time()
 
def numDigits(n):
    r"""
    return a sorted list of digits in a number
    """
    L = []
    while True:
        L.append(n % 10)
        n /= 10
        if n == 0: break
    L.sort()
    return L
 
i = 1
while True:
    a = numDigits(i)
    for j in range(1,7):
        b = numDigits(i * j)
        if a != b: break
        if j == 6:
            print i
            i = 0
    if i == 0: break
    i += 1
 
elapsed = time.time() - start
 
print "result found in %s seconds" % elapsed

When executed, we find the following result.

142857
result found in 1.92401385307 seconds

When I get into a slump (at work, while grading, etc.), one of my new-found activities to get myself rolling again is to find a Project Euler problem that I can solve before the cookie on the site expires. This one took about two minutes to form an algorithm for, another minute or so to code, and then 40 seconds to execute. I normally try to optimize these programs, but in this case I was just looking for a solution as fast as I could generate one.

Problem

145 is a curious number, as 1! + 4! + 5! = 1 + 24 + 120 = 145.
Find the sum of all numbers which are equal to the sum of the factorial of their digits.
Note: as 1! = 1 and 2! = 2 are not sums they are not included.

Path to a Solution

Let’s assume that an integer $d$ has $n$ digits and that $d$ is the sum of the factorials of those $n$ digits. How big can $n$ be? We quickly find that $n<8$ since the largest sum of factorials we can produce from an 8-digit number is itself a 7-digit number (set all digits equal to 9 and look at $8\times 9!$, which only has 7 digits). So, we know the numbers that we're looking at are at most 7-digit numbers. Since $9999999$ would give $7\times 9!=2,540,160$, we have the following implementation in Sage:

Sage Solution

import time
 
start = time.time()
 
# cache the factorials
F = [factorial(ZZ(i)) for i in range(0,10,1)]
 
S = 0
 
for i in range(3,2540161,1):
    # form sums of factorials
    D = sum([F[j] for j in ZZ(i).digits()])
    # is the sum of factorials equal to original number?
    if D == i: S += D
 
elapsed = time.time() - start
 
print "Result %s found in %s seconds." % (S, elapsed)

When executed, we get the correct result.

Result 40730 found in 41.8051738739 seconds.

In fact, there are only two numbers that have this property: 145 and 40585. In the end, the only thing I actually got from using Sage is the factorial function (and I also have to force the iterator values to be Sage integers).


The decimal number, $585=1001001001_2$ (binary), is palindromic in both bases. Find the sum of all numbers, less than one million, which are palindromic in base 10 and base 2. (Please note that the palindromic number, in either base, may not include leading zeros.)

Sage Solution

This is relatively straightforward in Sage.

import time
 
def is_palindromic(n):
    d = n.digits()
    if len(d) == 1: return True
    for i in range(floor(len(d)/2)):
        if d[i] != d[len(d)-1-i]: return False
    return True
 
def is_b_palindromic(n):
    n = ZZ(bin(n).split('0b')[1])
    return is_palindromic(n)
 
start = time.time()
 
s = 0
for n in range(1000000):
    if is_palindromic(ZZ(n)) and is_b_palindromic(ZZ(n)): s += n
 
elapsed = time.time() - start
 
print "result %s found in %s seconds" % (s, elapsed)

Executing that bit of code, we get the following result.

result 872187 found in 21.3442850113 seconds

A triangular pyramid is constructed using spherical balls so that each ball rests on exactly three balls of the next lower level.

Then, we calculate the number of paths leading from the apex to each position: A path starts at the apex and progresses downwards to any of the three spheres directly below the current position. Consequently, the number of paths to reach a certain position is the sum of the numbers immediately above it (depending on the position, there are up to three numbers above it).

The result is Pascal’s pyramid and the numbers at each level n are the coefficients of the trinomial expansion $(x + y + z)^n$. How many coefficients in the expansion of $(x + y + z)^{200000}$ are multiples of $10^{12}$?

Solution Using the Multinomial Theorem

The generalization of the binomial theorem is the multinomial theorem. It says that multinomials raised to exponents can be expanded using the formula
\[
(x_1+x_2+\cdots+x_m)^n=\sum_{{k_1+k_2+\cdots+k_m=n}\atop{0\le k_i\le n}}\left({n}\atop{k_1,k_2,\ldots,k_m}\right)\prod_{1\le t\le m}x_t^{k_t}
\]
where
\[
\left({n}\atop{k_1,k_2,\ldots,k_m}\right)=\frac{n!}{k_1!k_2!\cdots k_m!}.
\]
Of course, when m=2 this gives the binomial theorem. The sum is taken over all partitions $k_1+k_2+\cdots+k_m=n$ for integers $k_i$. If n=200000 abd m=3, then the terms in the expansion are given by
\[
\left({200000}\atop{k_1,k_2,k_3}\right)x_1^{k_1}x_2^{k_2}x_3^{k_3}=\frac{200000!}{k_1!k_2!k_3!}x_1^{k_1}x_2^{k_2}x_3^{k_3}
\]
where $k_1+k_2+k_3=200000$. It’s worth pointing out that each of the coefficients is an integer, and thus has a unique factorization into products of prime integers. Of course, there’s no way that we’re going to calculate these coefficients. We only need to know when they’re divisible by $10^{12}$. Thus, we only need to consider how many factors of 2 and 5 are involved.

First, we’ll create a function $p(n,d)$ that outputs how many factors of $d$ are included in $n!$. We have that
\[
p(n,d)=\left\lfloor\frac{n}{d}\right\rfloor+\left\lfloor\frac{n}{d^2}\right\rfloor+\left\lfloor\frac{n}{d^3}\right\rfloor+ \cdots+\left\lfloor\frac{n}{d^r}\right\rfloor,
\]
where $d^r$ is the highest power of $d$ dividing $n$. For instance, there are 199994 factors of 2 in 200000!. Since we’re wondering when our coefficients are divisible by $10^{12}=2^{12}5^{12}$, we’ll be using the values provided by $p(n,d)$ quite a bit for $d=2$ and $d=5$. We’ll store two lists:
\[
p2=[p(i,2)\text{ for }1\le i\le 200000]\quad\text{and}\quad p5=[p(i,5)\text{ for }1\le i\le 200000].
\]
For a given $k_1,k_2,k_3$, the corresponding coefficient is divisible by $10^{12}$ precisely when
\[
p2[k_1]+p2[k_2]+p2[k_3]<199983\ \text{and}\ p5[k_1]+p5[k_2]+p5[k_3]<49987.
\]
That is, this condition ensures that there are at least 12 more factors of 2 and 5 in the numerator of the fraction defining the coefficients.

Now, we know that $k_1+k_2+k_3=200000$, and we can exploit symmetry and avoid redundant computations if we assume $k_1\le k_2\le k_3$. Under this assumption, we always have
\[
k_1\le\left\lfloor\frac{200000}{3}\right\rfloor=66666.
\]
We know that $k_1+k_2+k_3=200000$ is impossible since 200000 isn't divisible by 3. It follows that we can only have (case 1) $k_1=k_2 < k_3$, or (case 2) $k_1 < k_2=k_3$, or (case 3) $k_1 < k_2 < k_3$.

In case 1, we iterate $0\le k_1\le 66666$, setting $k_2=k_1$ and $k_3=200000-k_1-k_2$. We check the condition, and when it is satisfied we record 3 new instances of coefficients (since we may permute the $k_i$ in 3 ways).

In case 2, we iterate $0\le k_1\le 66666$, and when $k_1$ is divisible by 2 we set $k_2=k_3=\frac{200000-k_1}{2}$. When the condition holds, we again record 3 new instance.

In case 3, we iterate $0\le k_1\le 66666$, and we iterate over $k_2=k_1+a$ where $1\le a < \left\lfloor\frac{200000-3k_1}{2}\right\rfloor$. Then $k_3=200000-k_1-k_2$. When the condition holds, we record 6 instances (since there are 6 permutations of 3 objects).

Cython Solution

I’ll provide two implementations, the first written in Cython inside Sage. Then, I’ll write a parallel solution in C.

%cython
 
import time
from libc.stdlib cimport malloc, free
 
head_time = time.time()
 
cdef unsigned long p(unsigned long k, unsigned long d):
    cdef unsigned long power = d
    cdef unsigned long exp = 0
    while power <= k:
        exp += k / power
        power *= d
    return exp
 
cdef unsigned long * p_list(unsigned long n, unsigned long d):
    cdef unsigned long i = 0
    cdef unsigned long * powers = <unsigned long *>malloc((n+1)*sizeof(unsigned long))
    while i <= n:
        powers[i] = p(i,d)
        i += 1
    return powers
 
 
run_time = time.time()
 
# form a list of number of times each n! is divisible by 2.
cdef unsigned long * p2 = p_list(200000,2)
 
# form a list of number of times each n! is divisible by 5.
cdef unsigned long * p5 = p_list(200000,5)
 
cdef unsigned long k1, k2, k3, a
cdef unsigned long long result = 0
 
k1 = 0
while k1 <= 66666:
    # case 1: k1 = k2 < k3
    k2 = k1
    k3 = 200000 - k1 - k2
    if 199982 >= (p2[k1]+p2[k2]+p2[k3]) and 49986 >= (p5[k1]+p5[k2]+p5[k3]):
        result += 3
    # case 2: k1 < k2 = k3
    if k1 % 2 == 0:
        k2 = (200000 - k1)/2
        k3 = k2
        if 199982 >= (p2[k1]+p2[k2]+p2[k3]) and 49986 >= (p5[k1]+p5[k2]+p5[k3]):
            result += 3
    # case 3: k1 < k2 < k3
    a = 1
    while 2*a < (200000 - 3*k1):
        k2 = k1 + a
        k3 = 200000 - k1 - k2
        if 199982 >= (p2[k1]+p2[k2]+p2[k3]) and 49986 >= (p5[k1]+p5[k2]+p5[k3]):
            result += 6
        a += 1
    k1 += 1
 
 
free(p2)
free(p5)
 
 
elapsed_run = round(time.time() - run_time, 5)
elapsed_head = round(time.time() - head_time, 5)
 
print "Result: %s" % result
print "Runtime: %s seconds (total time: %s seconds)" % (elapsed_run, elapsed_head)

When executed, we find the correct result relatively quickly.

Result: 479742450
Runtime: 14.62538 seconds (total time: 14.62543 seconds)

C with OpenMP Solution

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <omp.h>
 
/*****************************************************************************/
/* function to determine how many factors of 'd' are in 'k!'                 */
/*****************************************************************************/
unsigned long p(unsigned long k, unsigned long d) {
    unsigned long power = d;
    unsigned long exp = 0;
    while (power <= k) {
        exp += k/power;
        power *= d;
    }
    return exp;
}
 
/*****************************************************************************/
/* create a list [p(0,d),p(1,d),p(2,d), ... ,p(n,d)] and return pointer      */
/*****************************************************************************/
unsigned long * p_list(unsigned long n, unsigned long d) {
    unsigned long i;
    unsigned long * powers = malloc((n+1)*sizeof(unsigned long));
    for (i=0;i<=n;i++) powers[i] = p(i,d);
    return powers;
}
 
/*****************************************************************************/
/* main                                                                      */
/*****************************************************************************/
int main(int argc, char **argv) {
    unsigned long k1, k2, k3, a;
    unsigned long long result = 0;
 
    unsigned long * p2 = p_list(200000, 2);
    unsigned long * p5 = p_list(200000, 5);
 
 
    #pragma omp parallel for private(k1,k2,k3,a) reduction(+ : result)
    for (k1=0;k1<66667;k1++) {
        // case 1: k1 = k2 < k3
        k2 = k1;
        k3 = 200000 - k1 - k2;
        if (p2[k1]+p2[k2]+p2[k3]<199983 && p5[k1]+p5[k2]+p5[k3]<49987) {
            result += 3;
        } 
        // case 2: k1 < k2 = k3
        if (k1 % 2 == 0) {
            k2 = (200000 - k1)/2;
            k3 = k2;
            if (p2[k1]+p2[k2]+p2[k3]<199983 && p5[k1]+p5[k2]+p5[k3]<49987) {
                result += 3;
            }
        }
        // case 3: k1 < k2 < k3
        for (a=1;2*a<(200000-3*k1);a++) {
            k2 = k1 + a;
            k3 = 200000 - k1 - k2;
            if (p2[k1]+p2[k2]+p2[k3]<199983 && p5[k1]+p5[k2]+p5[k3]<49987) {
                result += 6;
            }
        }
    }
 
    free(p2);
    free(p5);
 
    printf("result: %lld\n", result);
 
    return 0;
}

This can be compiled and optimized using GCC as follows.

$ gcc -O3 -fopenmp -o problem-154-omp problem-154-omp.c

When executed on a 16-core machine, we get the following result.

$ time ./problem-154-omp 
result: 479742450
 
real    0m1.487s

This appears to be the fastest solution currently known, according to the forum of solutions on Project Euler. The CPUs on the 16-core machine are pretty weak compared to modern standards. When running on a single core on a new Intel Core i7, the result is returned in about 4.7 seconds.


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)):
    L[i].sort()
L.sort()
 
# 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:
    N.append([])
for p in P:
    s = p.digits()
    s.sort()
    if s in M:
        r = M.index(s)
        N[r].append(p)
 
# 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

Problem

Euler published the remarkable quadratic formula:

\[n^2+n+41\]

It turns out that the formula will produce 40 primes for the consecutive values $n=0$ to $39$. However, when $n=40$, $40^2+40+41=40(40+1)+41$ is divisible by 41, and certainly when $n=41$, $41^2+41+41$ is clearly divisible by 41.

Using computers, the incredible formula $n^2-79n+1601$ was discovered, which produces 80 primes for the consecutive values $n=0$ to $79$. The produce of the coefficients, $-79$ and $1601$ is $-126479$.

Considering quadratics of the form:

$n^2+an+b$ where $|a|

where $|n|$ is the modulus/absolute value of $n$, e.g., $|11|=11$ and $|-4|=4$.

Find the product of the coefficients $a$ and $b$, for the quadratic expression that produces the maximum number of primes for consecutive values of $n$, starting with $n=0$.

Sage Solution

Some observations:

  • Clearly, $b$ must be a prime or a negative prime, since $n=0$ must result in a prime.
  • We can refine a lower bound on $b$ as follows: If a given combination $(a,b)$ results in $m$ consecutive primes, then clearly $b>m$, since otherwise we would obtain a factor of $m$ in the given polynomial. Also, when $a$ is negative, we must have $b>-(n^2+an)$, since the prime values $n^2+an+b$ must be positive. Since we know that $n^2+n+41$ returns 40 primes, we know that any interesting values for $a$ and $b$ must then satisfy $b>-(40^2+40a)$.

Given these observations, I’ll write a routine in Sage.

import time
 
start = time.time()
 
P = prime_range(1000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
    for b in P:
        if b &lt; -1600-40*a or b &lt; c_max: continue
        c, n = 0, 0
        while is_prime(n**2 + a*n + b):
            c += 1
            n += 1
        if c &gt; c_max:
            a_max, b_max, c_max = a, b, c
 
prod = a_max * b_max
 
elapsed = time.time() - start
 
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

When executed, we obtain the following.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 8.53578400612 seconds

It’s entirely possible that we’re being a bit redundant with the integer values that we’re throwing at Sage’s “is_prime” function. We’ll cache a list of primes within the range that we’re considering. (I came up with the number 751000 below in the Python section.)

import time
 
start = time.time()
 
P = prime_range(751000)
L = [False] * 751000
for p in P: L[p] = True
 
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
    for b in range(1,1001):
        if L[b] is False: continue
        if b &lt; -1600-40*a or b &lt; c_max: continue
        c, n = 0, 0
        while L[n**2 + a*n + b] is True:
            c += 1
            n += 1
        if c &gt; c_max:
            a_max, b_max, c_max = a, b, c
 
prod = a_max * b_max
 
elapsed = time.time() - start
 
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

This indeed executes more quickly, as is seen in the output.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.39028286934 seconds

Python Solution

Of course, Sage’s prime/factorization functionality isn’t available in Python, and so we’re going to need to create a prime sieve and store primes in some form. One of the things I need to know is how large the values $n^2+an+b$ can be. If we were to assume that the maximum value of the length of a prime sequence is 500 (which is an overestimate as we now know), then we’d have the maximum value of the quadratic given as:

500**2 + 1000 * 500 + 1000
751000

Let’s create a prime sieve that will return a list of all primes below a given number.

def primes_xrange(stop):
    primes = [True] * stop
    primes[0], primes[1] = [False] * 2
    L = []
    for ind, val in enumerate(primes):
        if val is True:
            primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
            L.append(ind)
    return L

Let’s verify that it works.

primes_xrange(20)
[2, 3, 5, 7, 11, 13, 17, 19]

How long does this function take to form all of the primes under 751000?

time P = primes_xrange(751000)
Time: CPU 0.39 s, Wall: 0.39 s

Now, we’ll take the Sage routine from above and write it completely in Python. Here, I’m using only the sieve list portion of the primes_xrange function I just created. As above, I’m just checking to see if values in that list are True or False, corresponding to primes and composite numbers. This should take a bit longer to run than the Sage version above.

import time
 
def primes_xrange(stop):
    primes = [True] * stop
    primes[0], primes[1] = [False] * 2
    for ind, val in enumerate(primes):
        if val is True:
            primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
    return primes
 
start = time.time()
 
P = primes_xrange(751000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
    for b in range(1,1001):
        if P[b] is False: continue
        if b &lt; -1600-40*a or b &lt; c_max: continue
        c, n = 0, 0
        while P[n**2 + a*n + b] is True:
            c += 1
            n += 1
        if c &gt; c_max:
            a_max, b_max, c_max = a, b, c
 
prod = a_max * b_max
 
elapsed = time.time() - start
 
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

And it does run slightly slow compared to the Sage version, as expected.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.73107814789 seconds

Cython Solution

I’m going to take roughly the same approach as with the Python solution, but I’m going to malloc a C array instead of using the Python list. That should make things a bit more efficient. Before we do that, let’s see how much faster our routine runs if we just slap a “%cython” on the first line.

%cython
 
import time
 
def primes_xrange(stop):
    primes = [True] * stop
    primes[0], primes[1] = [False] * 2
    for ind, val in enumerate(primes):
        if val is True:
            primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
    return primes
 
start = time.time()
 
P = primes_xrange(751000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
    for b in range(1,1001):
        if P[b] is False: continue
        if b &lt; -1600-40*a or b &lt; c_max: continue
        c, n = 0, 0
        while P[n**2 + a*n + b] is True:
            c += 1
            n += 1
        if c &gt; c_max:
            a_max, b_max, c_max = a, b, c
 
prod = a_max * b_max
 
elapsed = time.time() - start
 
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

We find that this does improve the timing, and our Python code (which hasn’t been modified at all) now runs more quickly than the Sage code.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.02033305168 seconds

Rewriting to take advantage of Cython’s C datatypes, how much faster can we obtain the result?

%cython
 
import time
from libc.stdlib cimport malloc, free
 
start = time.time()
 
cdef long i, j
cdef bint *primes = malloc(751000 * sizeof(bint))
primes[0], primes[1] = 0, 0
i = 2
while i &lt; 751000:
    primes[i] = 1
    i += 1
i = 2
while i &lt; 751000:
    if primes[i] == 1:
        j = 2 * i
        while j &lt; 751000:
            primes[j] = 0
            j += i
    i += 1
 
cdef long long a, b, c, t
cdef long long pol, n
cdef long a_max, b_max, c_max, prod
a_max, b_max, c_max = 0, 0, 0
a = -1000
while a &lt; 1001:
    b = 2
    while b &lt; 1001:
        if primes[b] == 0:
            b += 1
            continue
        t = -1600-40*a
        if b &lt; t or b &lt; c_max:
            b += 1
            continue
        c, n = 0, 0
        pol = n**2 + a*n + b
        while primes[pol] == 1:
            c += 1
            n += 1
            pol = n**2 + a*n + b
        if c &gt; c_max:
            a_max, b_max, c_max = a, b, c
        b += 1
    a += 1
 
prod = a_max * b_max
 
free(primes)
 
elapsed = time.time() - start
 
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

Now, this should run much more quickly.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 0.0381400585175 seconds

Thus, the Sage version runs 1.245 times as fast as the Python version. The Cython version runs 45.387 times as fast as the Python version and 36.452 times as fast as the Sage version.


$n!$ means $n\times(n-1)\times\cdots\times 3\times 2\times 1$. For example, $10!=10\times 9\times\cdots\times 2\times 1=3628800$ and the sum of the digits in the number $10!$ is 27. Find the sum of the digits in $100!$.

Sage Solution

There is a single-line solution in Sage. I’ve added a bit more here just to determine the timing.

import time
 
start = time.time()
 
s = sum(factorial(100).digits())
 
elapsed = time.time() - start
 
print "Solution %s found after %s seconds" % (s,elapsed)

When executed, we find the following.

Solution 648 found after 0.000201940536499 seconds

Python Solution

Without the built-in functions of Sage, we can build a factorial function and use some string functionality in Python to achieve the same result.

import time
 
def factorial(n):
    if n == 0: return 1
    else: return n * factorial(n - 1)
 
def sum_digits(n):
    return sum([int(i) for i in str(n)])
 
start = time.time()
 
sumfac = sum_digits(factorial(100))
 
elapsed = time.time() - start
 
print("%s found in %s seconds") % (sumfac,elapsed)

When executed, the same result is returned.

648 found in 0.000312805175781 seconds

We could also compute the factorial using a reduced lambda function as in the following code.

#!/usr/bin/python
 
import time
 
 
def factorial(n):
    return reduce(lambda x,y: x*y, range(1,n+1,1))
 
def sum_digits(n):
    return sum([int(i) for i in str(n)])
 
start = time.time()
 
sumfac = sum_digits(factorial(100))
 
elapsed = time.time() - start
 
print("%s found in %s seconds") % (sumfac,elapsed)

Problem: $2^{15}=32768$ and the sum of its digits is $3+2+7+6+8=26$.
What is the sum of the digits of the number $2^{1000}$?

Sage Solution

Sage’s built-in Python and functions makes this easy.

import time
 
start = time.time()
a = 2^1000
s = sum(a.digits())
elapsed = time.time() - start
 
print "%s found in %s seconds" % (s,elapsed)

It executes pretty quickly too.

1366 found in 0.000343084335327 seconds

An Easy Python Solution

Python itself also makes this problem too easy, due to string functions.

import time
 
def pow2sum(exp):
    pow = list(str(2**1000))
    return sum([int(i) for i in pow])
 
start = time.time()
n = pow2sum(1000)
elapsed = (time.time() - start)
print "%s found in %s seconds" % (n,elapsed)

And it is fast:

1366 found in 0.000911951065063 seconds

Python Without Strings Attached

Let’s do our arithmetic without string functionality. Then, I’d note that $2^{1000}<10^{1000}$ and so we know that, at most, we're dealing with 1001 digits. So, we can create a routine to multiply 2 by itself 1000 times, maintaining the result at each step in a list instead of a single integer. (This is how you'd be forced to do things in C, where arbitrary length integers are pure fiction.) Since we're only multiplying by 2 at each iteration, we know that we'll either carry a zero or one to the next stage... which does make this routine a bit more simple than your typical multiply-by-list routine.

import time
 
def pow2sum(exp):
    L = [0] * exp # make a list exp long
    L[0] = 1
    for power in range(exp):
        carry, add = 0, 0
        for index in range(exp):
            prod = L[index] * 2 + carry
            if prod > 9:
                carry = 1
                prod = prod % 10
            else: carry = 0
            L[index] = prod
    return sum(L)
 
start = time.time()
n = pow2sum(1000)
elapsed = (time.time() - start)
print "%s found in %s seconds" % (n,elapsed)

It runs relatively quickly, although not as quickly as the string version runs.

1366 found in 0.361020803452 seconds

Cython Solutions

We can first trivially rewrite the string Python version.

%cython
 
import time
 
cdef pow2sum(unsigned int exp):
    cdef list pow = list(str(2**1000))
    return sum([int(i) for i in pow])
 
start = time.time()
n = pow2sum(1000)
elapsed = (time.time() - start)
print "%s found in %s seconds" % (n,elapsed)

When executed, we find that the string Cython code runs about 1.14 times as fast as the string Python code.

1366 found in 0.000799894332886 seconds

Of course, I’m more interested in seeing how much faster the arithmetic version runs.

%cython
 
import time
from libc.stdlib cimport malloc, free
 
cdef pow2sum(unsigned int exp):
    cdef int *L = <int *>malloc(exp * sizeof(int))
    cdef unsigned int power = 0,index = 0
    cdef unsigned int prod, carry, add
    while index < exp:
        L[index] = 0
        index += 1
    L[0] = 1
    while power < exp:
        carry, add, index = 0, 0, 0
        while index < exp:
            prod = L[index] * 2 + carry
            if prod > 9:
                carry = 1
                prod = prod % 10
            else: carry = 0
            L[index] = prod
            index += 1
        power += 1
    cdef int sum = 0
    index = 0
    while index < exp:
        sum += L[index]
        index += 1
    free(L)
    return sum
 
start = time.time()
n = pow2sum(1000)
elapsed = (time.time() - start)
print "%s found in %s seconds" % (n,elapsed)

When executed, we get the following result.

1366 found in 0.00858902931213 seconds