Euler-Project(I)

1. Multiples of 3 and 5

Problem Description

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

Solution

# below 1000 
# multiples of 3 and 5
# 3*(1+2+3+...+999/3) + 5*(1+2+3+...+999/5) - 15*(1+2+3+...+999/15)
# Note that 1+2+3+...+p = 1/2*p*(p+1)

count_result = 0

for item in range(3, 1000, 3):
    count_result += item

for item2 in range(5, 1000, 5):
    if item2 % 3 == 0:
        continue
    else:
        count_result += item2

print count_result

2. Even Fibonacci numbers

Problem Description

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

Solution

# 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
# odd even odd odd even odd odd even odd odd even
# Find the total number of the Fibonacci numbers below four million
origin1 = 1
origin2 = 2
count = 0
final_result = 2
origin_new = 0
while origin_new < 4000000:
    count += 1
    origin_new = origin1 + origin2
    if count % 3 == 0:
        final_result += origin_new
    origin1, origin2 = origin2, origin_new
print final_result

\1.png)

\2.png)

3. Largest prime factor

Problem Description

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

Solution1

# Largest prime factor of designated number 600851475143
# Note the key word: prime factor
import math
def isPrime(n):
    '''
    Those numbers who are some times of evens or some odd must be composite numbers
    '''
    if n == 2:
        return True
    elif n % 2 == 0:
        return False
    else:
        for item in range(3, int(math.sqrt(n))+1, 2):
            if n % item == 0:
                return False
    return True
def isPrime2(n):
    '''
    All numbers can be shown as such numbers:6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5
    In these numbers, besides 2 and 3, only 6n+1 and 6n+5 may be primes.
    For those who can be shown by 6n+1 and 6n+5 but aren't primes, they can be some times of 6n+1 or 6n+5
    '''
    if n == 1: return False
    if n == 2 or n == 3:
        return True
    if n % 6 != 1 and n % 6 != 5:
        return False
    for i in range(int(math.sqrt(n))+1):
        division = 5 + (i+1)*6 
        if (n // division > 1) or (n // (division+2) > 1):
            # division == 6n+5, division+2 == 6n+1
            return False
    return True
possible_result = []
final_number = 600851475143
item = int(math.sqrt(final_number))+1
while item >= 1:
    if final_number % item == 0:
        possible_result.append(item)
    item -= 2
for result in possible_result:
    if isPrime(result):
        print result
        break

Solution2

本题也可以使用厄拉多塞筛法寻找输入数字范围内所有的素数。因为素数的倍数一定不是素数，所以我们找到一个素数时可以将其倍数从所给范围内排除。这种方法称为素数筛。例如求100以内的素数：

n = 100
L = list(range(2, n))
ans = set()
while L:
    x = L.pop(0)
    ans.add(x)
    i = 2
    while i*x < n:
        if i*x in L:
            L.remove(i*x)
        i += 1
print(ans)
# ans = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}

基于以上思想，在找输入数字因数时可以将合数筛出，代码如下：

ans = []
n = 600851475143
iter_max = int(n ** 0.5)
for num in range(2,iter_max):
    if n%num == 0:
        ans.append(num)
        n/=num
        while n%num == 0:
            n/=num # 保证n已被num除尽，此时n不会再有num*i的因数
print(ans)
ans = [71, 839, 1471, 6857]

素数判别法

根号判别法

设输入的数字为n，则可以通过遍历方法暴力搜索其素因子，若出现非1及其本身的素因子，则可断定该数字为素数。常见的遍历范围为1—n-1，其实将遍历范围调节至1—$\sqrt{n}$亦可实现素数判定的目的。

奇偶判别法

对于所有可能成为数字x素因子的n-1个数字，偶数中除了2均不是质数，且奇数的因数没有偶数，因此可以继续优化。首先将n与2进行比较，其次判断2是否为n的素因子，最后从3开始以2为增幅逐次判断数字n是否包含奇数因子。

def isPrime(n):
    '''
    Those numbers who are some times of evens or some odd must be composite numbers
    '''
    if n == 2:
        return True
    elif n % 2 == 0:
        return False
    else:
        for item in range(3, int(math.sqrt(n))+1, 2):
            if n % item == 0:
                return False
    return True

6n系判别

所有数字均可表示为6n，6n+1，6n+2，6n+3，6n+4，6n+5的形式，除2和3以外，所有的素数都可以表示为6n+1和6n+5的形式，如果输入数字x是6n+1和6n+5的整数倍，则必为合数。

def isPrime2(n):
    '''
    All numbers can be shown as such numbers:6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5
    In these numbers, besides 2 and 3, only 6n+1 and 6n+5 may be primes.
    For those who can be shown by 6n+1 and 6n+5 but aren't primes, they can be some times of 6n+1 or 6n+5
    '''
    if n == 1: return False
    if n == 2 or n == 3:
        return True
    if n % 6 != 1 and n % 6 != 5:
        return False
    for i in range(int(math.sqrt(n))+1):
        division = 5 + (i+1)*6 
        if (n // division > 1) or (n // (division+2) > 1):
            # division == 6n+5, division+2 == 6n+1
            return False
    return True

Miller-Rabin素数判别法

该素数判别方法应用费马小定理对素数进行概率判定，若输入数字N通过t次测试，则N不是素数的概率仅为$(1/4)^{t}$，随着通过测试次数的增加，N是素数的概率无穷逼近于1。在实际运用中，可首先用300—500个小素数对N进行测试，以提高测试通过的概率与算法的速度。

具体步骤如下：

• 计算奇数M，使得N=$2^{r}*M+1$；
• 选择随机数A<N；
• 对于任意i<r，若$A^{(2^{i}*M)}mod N=N-1$，则N通过随机数A的测试；
• 或者若$A^{M}modN=1$，则N通过随机数A的测试；
• 改变随机数A的值对N进行多次测试（一般为5—10次，较高需求的情况下可进行20—30次），若全部通过则判定N为素数。

相关代码如下：

import random
def fast_power(base, power, n):
    result = 1
    tmp = base
    while power > 0:
        if power&1 == 1:
            result = (result * tmp) % n
        tmp = (tmp * tmp) % n 
        power = power>>1
    return result
def Miller_Rabin(n, iter_num):
    # 2 is prime
    if n == 2:
        return True
    # if n is even or less than 2, then n is not a prime
    if n&1 == 0 or n<2:
        return False
    # n-1 = (2^s)m
    # get random odd m
    m,s = n - 1,0
    while m&1==0:
        m = m>>1
        s += 1
    # M-R 
    for _ in range(iter_num):
        '''
        key algorithm
        '''
        b = fast_power(random.randint(2,n-1), m, n)
        if b==1 or b== n-1:
            continue
        for __ in range(s-1):
            b = fast_power(b, 2, n)
            if b == n-1:
                break
        else:
            return False
    return True
if __name__ == "__main__":
    # example
    print(Miller_Rabin(49139, 10))
    print(Miller_Rabin(561, 10))

素数筛查

除了上述的高级试除法外，还可以使用素数筛查的方法。常见的素数筛查包括埃拉托斯特尼筛法和欧拉筛法。

埃氏筛法由希腊数学家埃拉托斯特尼提出，用以简单鉴定素数，方法如下：要获取自然数n（上界）内的全部素数，必须剔除所有小于等于sqrt(n)的素数的倍数，经过此种筛查后，剩下的就是素数。

欧拉筛法是埃氏筛法的改进。采用欧拉筛法进行筛选过程中，对于含多个因子的数字需要进行多次筛选，耗费部分运行时间。例如，对于合数20，可分解为2*10，4*5，故至少需要筛选两次。欧拉筛过程中引入语句if i%prime[j] == 0: break，保证每个合数只被这个合数最小的质因子筛除，而且不出现重复筛除。

代码实现可参考：素数筛法详解（欧拉筛&埃氏筛）。

4. Largest palindrome product

Problem Description

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.

Find the largest palindrome made from the product of two 3-digit numbers.

Solution

'''
    Suppose that P = a*b = 100000*x+10000*y+1000*z+100*z+10*y+x = 11*(9091*x+910*y+100*z)
    The range of a and b are both 100 to 999
'''
def reverse(input_number):
    reversed = 0
    while input_number > 0:
        reversed = reversed*10 + input_number % 10
        input_number /= 10
    return reversed

def isPalindrome(input_number):
    return input_number == reverse(input_number)

largest_number = 0
a = 999
while a >= 100:
    if a % 11 == 0:
        # Then b doesn't need to contain prime factor 11
        b = 999
        b_down = 1
    else:
        # b needs to contain prime factor 11
        b = 990
        b_down = 11
    while b >= 100:
        if a*b <= largest_number:
            break
        if isPalindrome(a*b):
            largest_number = a*b
        b = b - b_down
    a = a - 1
print largest_number

5. Smallest multiple

Problem Description

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

Solution

'''
    The factors are 20 numbers from 1 to 20
'''
import math
def isPrime(n):
    '''
    All numbers can be shown as such numbers:6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5
    In these numbers, besides 2 and 3, only 6n+1 and 6n+5 may be primes.
    For those who can be shown by 6n+1 and 6n+5 but aren't primes, they can be some times of 6n+1 or 6n+5
    '''
    if n == 1: return False
    if n == 2 or n == 3:
        return True
    if n % 6 != 1 and n % 6 != 5:
        return False
    for i in range(int(math.sqrt(n))+1):
        division = 5 + (i+1)*6 
        if (n // division > 1) or (n // (division+2) > 1):
            # division == 6n+5, division+2 == 6n+1
            return False
    return True

factors_list = []
original_list = [x+1 for x in range(20)]
print original_list
# Get all primes
for item in original_list:
    if isPrime(item):
        factors_list.append(item)

k = 20
i = 0
N = 1
edge_number = math.sqrt(k)
while True:
    if i == len(factors_list):
        break
    if factors_list[i] <= edge_number:
        N = N*pow(factors_list[i], math.floor(math.log(k)/math.log(factors_list[i])))
        i += 1
        print N
    else:
        N = N*factors_list[i]
        i += 1
        print N
print N

分析题目，首先以1—10之间的数字为例，我们进行以下操作：

• Step1：寻找2—10内所有素数，即2、3、5、7，则剩余合数均可以素数的幂乘积的形式表示；待求的最小倍数可以表示为$min_multiply = 2^a3^b5^c*7^d$。
• Step2：分解剩余合数，均表示为min_multiply所示形式。
• Step3：取最大指数作为待求参数a，b，c，d的值。

现令N为可被2—k间所有数字整除的最小数字，分析求解N的过程。参考上述k=10时的分析，首先确定所有小于k的素数，存入列表P。随后确定列表中每个元素的次数，令$P[i]^{a[i]} = k$，两侧同时进行对数运算并向下取整，得$a[i] = floor(log(k) / log(P[i]))$。值得注意的是，当$P[i]^{2} > k$时，a[i]==1，故只需要计算满足$P[i] <= \sqrt(k)$对应素数的次数a[i]。最终$N=P[0]^{a[0]}P[1]^{a[1]}P[2]^{a[2]}*……..$

6. Sum square difference

Problem Description

The sum of the squares of the first ten natural numbers is,

$$1^2 + 2^2 + ... + 10^2 = 385$$

The square of the sum of the first ten natural numbers is,

$$(1 + 2 + ... + 10)^2 = 55^2 = 3025$$

Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is $3025 - 385 = 2640$.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

Solution

'''
    sum1 = 1+2+3+..+n = n*(n+1)/2
    sum2 = 1**2 + 2**2 + 3**2 +...+ n**2 = n*(n+1)*(2*n+1)/6
'''
sum1 = 100*101/2
sum2 = 100*101*201/6
print sum2-sum1

7. 10001st prime

Problem Description

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10 001st prime number?

Solution1

题目只给出素数的个数，最简单的思路就是循环计数与素数判断相结合。素数判断过程中，现补充事实如下：

• 1不是素数
• 除2以外的所有素数均为奇数
• 所有比3大的素数均可以写成6k+/-1的形式，k为整数
• 如果我们无法找到小于等于sqrt(n)并能够整除n的数字f，则可以判定n为素数

# Method1: Round and prime judgement
import math
def isPrime(input_number):
    if input_number == 1: return False
    elif input_number < 4: 
        return True 
    elif input_number % 2 == 0:
        return False
    elif input_number < 9:
        return True
    elif input_number % 3 == 0:
        return False
    else:
        limit = math.floor(math.sqrt(input_number))
        f = 5
        while f <= limit:
            # Which means 6n-1
            if input_number % f == 0: return False
            # We must plus 2 to get 6n+1
            if input_number % (f+2) == 0:
                return False
            f += 6
    return True
final_number = 3
count = 2
while count < 10001:
    final_number += 2
    if isPrime(final_number):
        count += 1
    else:
        continue
print final_number

Solution2

本题也可以采用埃拉托斯特尼筛法进行求解，求解的关键是确定第10001个素数的上界。为确保一定会出现第10001个素数，取较大的上界为1000000。

# Method2: a sieve of Eratosthenes
def findPrime(L):
    i = 0
    while L[i]**2 < L[-1]:
        # Eratosthenes algorithm
        count = 0
        for j in range(i+1, len(L)):
            if L[j] % L[i] == 0:
                L[j] = 0
                count += 1
        L.sort()
        L = L[count:]
        i += 1
    return L
# L = (2,3,4,...,1000000)
L = range(2, 1000000)
primes_list = findPrime(L)
print(primes_list[10000])

8. Largest product in a series

Problem Description

The four adjacent digits in the 1000-digit number that have the greatest product are 9 × 9 × 8 × 9 = 5832.

73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450

Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?

Solution

'''
    Find the thirteen adjacent digits in the 1000-digits number
'''
input_number = given_string
# Transform the 1000-digits number to number list in order to get each number seperately
input_number_list = [eval(input_number[i]) for i in range(len(input_number))]
maximum_number = 0
count = 0
while count <= len(input_number) - 13:
    # Combination of function reduce() and lambda function
    result_number = reduce(lambda x, y: x*y, input_number_list[count:count+13])
    if result_number > maximum_number:
        maximum_number = result_number
    count += 1
print maximum_number

9. Special Pythagorean triplet

Problem Description

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

$$a^{2} + b^{2} = c^{2}$$

For example, $3^{2} + 4^{2} = 9 + 16 = 25 = 5^{2}$.

There exists exactly one Pythagorean triplet for which a + b + c = 1000.
Find the product abc.

Solution1

可以采用放缩的方法，以题目所给的1000为例，设$ a=m1k，b=m2k，c=m3*k$，这里m1，m2，m3为小于50的勾股数。故只需要寻找满足1000 % (m1+m2+m3) == 0的勾股数并相应扩大k倍，使其满足a+b+c == 1000即可。

'''
    1. a+b+c = 1000
    2. Suppose a<b<c,then we get a^2 + b^2 = c^2
    3. Find a,b,c and calculate abc
'''
# First find all the Pythagorean triplet number under 50
for i in range(1, 50):
    for j in range(1, 50):
        for k in range(1, 50):
            if i**2+j**2 == k**2 and 1000%(i+j+k) == 0:
                print i*j*k*(1000/(i+j+k))**3
                break

Solution2

# Find the only Pythagorean triplet(a, b, c), for which a+b+c = 1000
# a^2+b^2=(s-a-b)^2, cause a<b<c, then a<=(s-3)/3 and b <(s-a)/2

def getPythagorean(s):
    for a in range(3, (s-3)//3):
        for b in range(a+1, (s-1-a)//2):
            c = s-a-b
            if c*c == a*a + b*b:
                print (a,b,c)
                return a*b*c

final_multiply = getPythagorean(1000)
print final_multiply

Solution3

下面通过勾股数的一些性质简化解法二中的代码，提高程序运行效率。

如果勾股数(a, b, c)满足gcd(a, b, c) = 1，定义该类勾股数具有素数性质。同时，任意组勾股数(a, b, c)可表示为：

$$a = m^{2}-n^{2}，b = 2*m*n，c = m^{2}+n^{2}，m>n>0 (9.1)$$

由于勾股数可通过倍乘进行放缩，故可对(a, b, c)进行运算，得到以下结果：

$$a = (m^{2}-n^{2})*d，b = 2*m*n*d，c = (m^{2}+n^{2})*d, m>n>0 (9.2)$$

使用以上性质，我们可得：

$$a+b+c = 2*m*(m+n)*d (9.3)$$

所以想要找到勾股数组合(a, b, c)满足a+b+c = s，我们需要在1—s/2之间寻找除数m，并寻找s/2m的奇除数k（此处k=m+n，k满足m < k < 2m，并且m，k互素）。随后令n = k - m，d = s/2mk，将该结果插入式(9.2)中，代码如下：

import math
def euclid_gcd(a, b):
    if a < b:
        a, b = b, a
    while b != 0:
        a, b = b, a%b
    return a

def getPythagorean(s2, mlimit):
    for m in range(2, mlimit):
        if s2 % m == 0:
            sm = s2 // m
            while sm % 2 == 0:
                sm = sm // 2
            if m % 2 == 1:
                k = m+2
            else:
                k = m+1
            while k < 2*m and k <= sm:
                if sm % k == 0 and euclid_gcd(k, m) == 1:
                    d = s2 // (k*m)
                    n = k - m
                    a = d*(m*m-n*n)
                    b = 2*d*m*n
                    c = d*(m*m+n*n)
                    print (a, b, c)
                    return a*b*c
                k += 2
s = 1000
s2 = s // 2
mlimit = int(math.ceil(math.sqrt(s2))) - 1
final_result = getPythagorean(s2, mlimit)
print final_result

10. Sum all primes below N million

Problem Description

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

Solution1

# Sum after the sieve of Eratosthenes
def findPrime(L):
    i = 0
    while L[i]**2 < L[-1]:
        # Eratosthenes algorithm
        count = 0
        for j in range(i+1, len(L)):
            if L[j] % L[i] == 0:
                L[j] = 0
                count += 1
        L.sort()
        L = L[count:]
        i += 1
    return L
L = range(2, 2000000)
primes_list = findPrime(L)
print reduce(lambda x,y:x+y, primes_list)

Solution2

除2以外的所有偶数均为合数，故只需要对所给范围2—N内的奇数进行素性判断。我们建立下标i与奇数2*i+1的对应关系。首先建立布尔类型数组，其长度为(N-1)/2，初始值均为False，表示全为素数。设p=2*i+1，则$p^{2}=4i^{2}+4i+1$，其对应的数组下标为2*i+1；设m=k*p，则m+2*p对应的下标为j+p。

同样参考埃氏筛的思路，下一步需要找外部循环的下标范围、内部循环的初始值以及步进参数。结合以上分析进行讨论。由$2i_max+1<= \sqrt(N)$可得i的范围为$(\lfloor N\rfloor -1) / 2$，该范围是外部循环的下标范围。对于奇数p，所有小于p^2的数字中若为合数则必可被小于p的素数整除，故内部循环的初始值设定为p^2，对应下标为2\i*(i+1)。又$p^{2}+p = 4i^{2}+4i+1+2i+1 = 4i^{2}+6i+2 = 2(2i^{2}+3i+1)$，可知若以p为步进，则得到的结果必为合数，无需进行判断。同时$p^{2}+2p = 2(2i+1)+1 = 4i^{2}+8*i+3 $，可以看出该结果必为合数，同时该结果为p的倍数，将数组对应下标位置标记为True。代码如下：

# Some imporvement of the sieve of Eratosthenes
# Only consider the odd numbers, do some index-arithmetics
import math
def findPrimes(sievebound, sieve_set, input_number):
    crosslimit = int(math.floor(math.sqrt(input_number)-1 // 2))
    for i in range(1, crosslimit):
        if not sieve_set[i]:
            for j in range(2*i*(i+1), sievebound, 2*i+1):
                # 2*j+1 is compositive
                sieve_set[j] = True
    return sieve_set
input_number = 2000000
sievebound = int(math.floor((input_number-1) // 2))
sieve_set = [False]*input_number
sieve_set_after_deal = findPrimes(sievebound, sieve_set, input_number)
sum_result = 2
for k in range(1, sievebound):
    if not sieve_set[k]:
        sum_result += 2*k+1
print sum_result

以题中所给的2000000为例，解法一运行速度为8s左右，解法二运行速度大大提升，只需要0.417s。

\3.png)