Performance of modular square root

Question

Here's Project Euler problem 451:

Consider the number \$15\$. There are eight positive numbers less than \$15\$ which are coprime to \$15\$: \begin{align*} 1, 2, 4, 7, 8, 11, 13, 14. \end{align*} The modular inverses of these numbers modulo \$15\$ are: \begin{align*} 1, 8, 4, 13, 2, 11, 7, 14 \end{align*} because

\begin{align*} 1 \times 1 \mod 15 &= 1 \\ 2 \times 8 = 16 \mod 15 &= 1 \\ 4 \times 4=16 \mod 15 &= 1 \\ 7 \times 13=91 \mod 15 &= 1 \\ 11 \times 11=121 \mod 15 &= 1 \\ 14 \times 14=196 \mod 15 &= 1 \end{align*}

Let \$I(n)\$ be the largest positive number \$m\$ smaller than \$n − 1\$ such that the modular inverse of \$m\$ modulo \$n\$ equals \$m\$ itself. So \$I(15) = 11\$. Also \$I(100) = 51\$ and \$I(7) = 1\$.

Find \$\sum I(n)\$ for \$3 \leq n \leq 2 \times 10^7\$.

I realised that because they are symmetrical, we can search for the smallest ones (of course skipping the trivial case \$1^2 = 1 \mod n\$). I also realised that if \$n\$ is divisible by a prime, I don't need to test multiples of that prime. I currently devised a implementation for the first four primes, as I think they are the ones that have the most impact. Anyways, here's my code. However, I must warn you, that in trying to optimise for speed, this one is very memory inefficient (needs around 2.5 GB).

upto=2*10**7+1

a = [True] * upto
p = []

for n in range(2,upto):
    if a[n]:
        p.append(n)
        for k in range(2,(upto+n-1)//n):
            a[k*n] = False
p=set(p)


su=0
squ=list(map(lambda x: x*x, range(upto)))

print('primes and squares ready')

s2 = set(range(0,upto,2))# for testing divisibility
s3 = set(range(0,upto,3))
s5 = set(range(0,upto,5))
s7 = set(range(0,upto,7))

r=[]
r.append([])
r.append([])
r[0].append([])
r[0].append([])
r[1].append([])
r[1].append([])
r[0][0].append([])
r[0][0].append([])
r[0][1].append([])
r[0][1].append([])
r[1][0].append([])
r[1][0].append([])
r[1][1].append([])
r[1][1].append([])
r[0][0][0].append(range(2,upto))# for iterating
r[1][0][0].append(filter(lambda x: x not in s2, range(upto)))
r[0][1][0].append(filter(lambda x: x not in s3, range(upto)))
r[0][0][1].append(filter(lambda x: x not in s5, range(upto)))
r[0][0][0].append(filter(lambda x: x not in s7, range(upto)))

print('list pre generated')

r[1][1][0].append(filter(lambda x: x  in r[1][0][0][0], r[0][1][0][0]))
r[1][1][1].append(filter(lambda x: x  in r[0][0][1][0], r[1][1][0][0]))
r[0][1][1].append(filter(lambda x: x  in r[0][0][1][0], r[0][1][0][0]))
r[1][0][1].append(filter(lambda x: x  in r[0][0][1][0], r[1][0][0][0]))
r[1][0][1].append(filter(lambda x: x  in r[1][0][1][0], r[0][0][0][1]))
r[1][0][0].append(filter(lambda x: x  in r[1][0][0][0], r[0][0][0][1]))
r[1][1][0].append(filter(lambda x: x  in r[1][1][0][0], r[0][0][0][1]))
r[1][1][1].append(filter(lambda x: x  in r[0][0][0][1], r[1][1][1][0]))
r[0][1][1].append(filter(lambda x: x  in r[0][0][0][1], r[0][1][1][0]))
r[0][0][1].append(filter(lambda x: x  in r[0][0][1][0], r[0][0][0][1]))
r[0][1][0].append(filter(lambda x: x  in r[0][1][0][0], r[0][0][0][1]))

print('set up complete')

for n in range(3,upto):

    if n in p:
        su+=1
        print('testing', n)
    else:

        if n in s2:
            d2=1
        else:
            d2=0

        if n in s3:
            d3=1
        else:
            d3=0

        if n in s5:
            d5=1
        else:
            d5=0

        if n in s7:
            d7=1
        else:
            d7=0

        for x in r[d2][d3][d5][d7]:
            if squ[x]%n==1:
                su+=n-x
                break

print(su)

I have tried optimising the speed of division testing using sets, and pre-generated all the generators or lists (I wasn't sure about the performance in this case, so I tried both). I must say the code runs fast at first, checking the first 100,000 moduli in under a minute, however close to a million is slow substantially. Although cProfile doesn't time built-in functions I believe the majority of time is spent in the for x in r[d2][d3][d5][d7]: loop.

I'm not sure if the algorithm I'm using is wrong and faster ones exist, or if my code is simply inefficient (Project Euler claims the problem can be computed in under a minute, not sure if this is the case with Python though).

As this is a project to benefit my understanding of programming, and algorithms, please don't post full solutions but rather guidance and pointers.

Answer 1

I'm not really sure I understand your code; it's not really organized in a way that aids readability. Use functions, so that I don't need to spend a minute or more reading the first few lines of the program to understand that they should be

def computePrimesUpTo(n):
    sieve = [True] * n
    p = []
    for i in range(2, n):  # In Python2, use xrange instead
        if sieve[i]:
            p.append(i)
            for multiple in range(2 * i, n, i):
                sieve[multiple] = False
    return set(p)
primes = computePrimesUpTo(upTo)

Now at a glance, I can see that if n in primes is checking whether n is a prime number, and I can see that you're using the Sieve of Eratosthenes to compute them. (Note: p is not the worst offender here; in fact, it's the easiest to figure out. It's taken me several reads through your code to figure out what r is supposed to be.)

However, it looks like you're precomputing x*x for x in [0, 2*10^7], as well as lists of all multiples of 2, 3, 5, and 7 in the same range. This belies a fundamental misunderstanding of what is expensive in modern processors. x will be in a register, or in the worst case in the L1 cache, and squaring it will require a single instruction (less than a nanosecond). squ[x] will frequently not be in any of your CPU's memory caches; as a rule of thumb a cache miss costs 100 ns (And that's if it hasn't been swapped out of main memory onto disk! See this list of numbers every programmer should know). Not only that, but by reading from them, you'll push something else out of your cache, slowing down everything else. squ, s2, s3, s5, and s7 are likely slowing your program down by a noticeable amount.

I am less sure about your rs; they may help or hurt. However, the computation of r[1][1][0], etc. is not as fast as it might be; the in operation of python lists is O(n). Try implementing this function:

def list_intersection(a, b):
    """Computes the intersection of two sorted iterables.

    Returns a list; requires O(len(a) + len(b)) time.
    """
    pass

---

OK, now let's talk algorithm.

First of all, when computing l(n) there's no need to consider numbers less than or equal to sqrt(n). To take advantage of that, you might reverse the r[w][x][y][z] lists and iterate over them backward, with a check for each n that r[d2][d3][d5][d7][-1] > sqrt(n). This will give you a small benefit.

I think you had the right insight: If you're computing the coprimes separately for each n,

for c in [c for c in range(2,n) if is_coprime(c, n)]: if (c*c)%n == 1: break

is much slower than

for c in range(2,n): if (c*c)%n == 1: break

Filtering down range(2,n) is useful, but there's diminishing returns: filtering out factors of p gets you at most 1/p speed-up while computing l(n) with an additional cost that's O(upTo). Still, any test you do in your for x in r[...] loop will only slow you down (since if (x*x)%n == 1 is essentially as fast as you can get), so the only speedup you can gain is by reducing the size of the r lists.

Answer 2

You should find an algorithm finding the root square of 1 module n.

x^2 = 1 mod n

To do that:

• Factorize n in prime factor. n = p1^e1 * p2^e2 * ... * pk^ek
• Find the root for all prime xi^2 = 1 mod pi^ei. easy, it is 1 and -1 (or pi^ei - 1)
• Then use the Chinese theorem to find the solution for x (x = x1 mod x1^e1, x = ..., x = xk mod pk^ek (xi can take 2 value so you will have k^2 answer)

you just need to take the solution which meet your condition

EDIT:

I made a mistake about the power of prime root square result:

• 1 and -1 are the root square of xi^2 = 1 mod pi^ei but there may be other when e1 != 1
• for example: n= 16 we have x^2 = 1 mod 2^4 where 9 is also a solution too (9*9 = 81 = 1)

Answer 3

I think you are going in the wrong direction.

1) Instead of using your lengthy program block to initialize r you can use

r=[[[None,None],[None,None]],[[None,None],[None,None]]]
r[0][0][0]=range(2,upto)
...

of course it is also possible to use

r=[[[[],[]],[[],[]]],[[[],[]],[[],[]]]]   

this looks like more funny but is rather confusing so I prefer the first solution

2) what youd did

for n in range(2,upto):
    if a[n]:
        for k in range(2,(upto+n-1)//n):
            a[k*n] = False

this tooks 40 seconds on my notebook

what you want

for n in range(2,upto):
    if a[n]:
        for k in range(2*n,upto,n):
            a[k] = False

this tooks 30 seconds on my notebook. That difference is not essential but your loop is more complex and slower.

3) It does not make any sense to store the square numbers instead of calculating them. It may make sense to store the calculation of time consuming calculations, but not the result of a simple fast operation like squaring.

The same is true for % and //.

I did not get the idea of your algorithm and don't want to decipher the rest of your program.

---

So you should rewrite the program without using the lists s2,s3,s5,s7 and r. Also you should be aware that if you use functional programming there are a lot of readers that are not familiar with and that therefore won't review your code.