Problem #133 says:

A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k; for example, R(6) = 111111.

Let us consider repunits of the form \(R(10^{n})\).

Although R(10), R(100), or R(1000) are not divisible by 17, R(10000) is divisible by 17. Yet there is no value of n for which \(R(10^{n})\) will divide by 19. In fact, it is remarkable that 11, 17, 41, and 73 are the only four primes below one-hundred that can be a factor of R(10n).

Find the sum of all the primes below one-hundred thousand that will never be a factor of R(10n).

This solution works based on the Multiplicative Order of 10 mod p; where p is a prime number (for this problem a prime under 100,000). For the number to never be a divisor of a RepUnit in the form \(R(10^{n})\) then the above multiplicative order has to be something that can not be written in the form \(2^{i} * 5^{i}\). This method expands to any repunit no matter how big the value of n in \(R(10^{n})\) is.

The other method is to calculate a sufficiently large power of 10 such as \(10^{20}\) and use that in pythons pow(10,\(10^{20}\),p) which is faster but if you don’t choose a large enough power of ten you’ll get bad results after a certain point.

The second method is faster but I use the first in the below solution.

#!/usr/bin/python
from math import sqrt
#prime generator code aquired from
#http://code.activestate.com/recipes/366178-a-fast-prime-number-list-generator/
def Primes(n):
	if n==2: return [2]
	elif n<2: return []
	s=range(3,n+1,2)
	mroot = n ** 0.5
	half=(n+1)/2-1
	i=0
	m=3
	while m <= mroot:
		if s[i]:
			j=(m*m-3)/2
			s[j]=0
			while j<half:
				s[j]=0
				j+=m
		i=i+1
		m=2*i+3
	return [2]+[x for x in s if x]

primes = Primes(100000)

def Factor(n):
    global primes
    pmax = int(sqrt(n))+1
    factors = []
    for prime in primes:
        if prime > pmax:
            break
        if (n % prime ==0):
            while (n% prime ==0):
                factors = factors + [prime]
                n = n/prime
    if (n != 1):
        factors += [n]
    return factors

#n is the number to check
#form is the list of factors
#only true if n contains only factors found in form
def OfForm(n,form):
    factors = Factor(n)
    if (factors == []):
        return False
    for factor in factors:
        if (factor not in form):
            return False
    return True

#a^k = 1 (mod n).
def MultiplicativeOrder (a,n):
    k = 1
    if (a%n ==0):
        return 0
    r = pow(a,k,n)
    while (r != 1):
        k= k+1
        r = pow(a,k,n)
    return k

primelist = Primes(100000)
summation = 0

for p in primelist:
    if not OfForm(MultiplicativeOrder (10,p),[2,5]):
        print p
        summation += p

print summation