Is there any more optimal way to solve this idempotent equation (modulo n ring)?

Question

This is one of the problems on Project Euler:

If we calculate a^2 mod 6 for 0 <= a <= 5 we get: 0, 1, 4, 3, 4, 1.

The largest value of "a" such that a^2 mod 6 = a is 4. Let's call M(n) the largest value of a < n such that a^2 mod n = a. So M(6) = 4.

Find M(n) for 1 <=n <=10^7.

So far, this is what I have:

import time
start = time.time()
from math import sqrt

squares=[]
for numba in xrange(0,10000001/2+2):
    squares.append(numba*numba)
def primes1(n):
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(3,int(sqrt(n))+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]


tot=0

gor = primes1(10000001)

def factor1(n):
    '''Returns whether a number has more than 1 prime factor'''
    boo = False


    '''if n in gor:
        return True'''


    for e in xrange(0,len(gor)):



        z=gor[e]    

        if n%z==0:

            if boo:
                return False
            boo = True
        elif z*2>n:  
            break
    return True

for n in xrange(2,10000001):

    if factor1(n):

        tot+=1
    else:
        for a in xrange(int(sqrt(n))+1,n/2+1):
            if squares[a]%n==a:

                tot+=n+1-a
                break

print tot
print time.time()-start

I've tried this code for smaller cases and it works perfectly; however, it is way too slow to do 10^7 cases.

Currently, for n being less than 20000, it runs in about 8 seconds. When n is less than 90000, it runs in about 150 seconds.

As far as I can tell, for n is less than 10^7, it will run for many hours if not days.

I'm already using the sieve to generate prime numbers so that part is as fast as it can be, is there anything I can do to speed up the rest of the code?

I've already tried using different compiler like psyco, pypy, and shedskin. Psyco provides a minimal increase, shedskin speeds it up about 7 times but creates errors when large numbers occur, pypy speeds it up the most (about 20-30x the speed). But even then, it's still not fast enough for the amount of cases it has to go through.

Edit:

I added

squares=[]
for numba in xrange(0,10000001/2+2):
    squares.append(numba*numba)

This pre-generates all the squares of a before-hand so that I don't have to keep generating the same ones over and over again. Program became slightly faster but still not enough

Here's the original Project Euler problem: http://projecteuler.net/problem=407 — Ilia Sucholutsky, May 06 '13 at 23:05
Sure there's a better way. Hint: Can you reduce the number of candidates for `a` by using the prime factorisation of `n`? — Daniel Fischer, May 06 '13 at 23:14
Might be better asked on http://codereview.stackexchange.com/ — kindall, May 06 '13 at 23:14
I can reduce the candidates for a by iterating using values of prime factors of n, however, when I tested that, it ended up being about the same speed because of the additionoal operations I had to perform. — Ilia Sucholutsky, May 07 '13 at 22:20
Look at the equation and think. And solve it yourself, it is more satisfying than cheating. — starblue, May 09 '13 at 16:54

score 0 · Answer 1 · answered May 08 '13 at 13:29

This might depend on the size of N because of memory usage, but in smaller tests I found something of an improvement by precalculating the factor counts. So something like this:

factors = [0]*N
for z in gor:
  for n in xrange(1,N):
    m = z*n
    if m >= N: break
    factors[m] += 1

where N is 10000001, or whatever counter you're using.

Then instead of if factor1(n) you do if factors[n] < 2.

Is there any more optimal way to solve this idempotent equation (modulo n ring)?

1 Answers1