Euler-Project(III)

21. Amicable numbers

Problem Description

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

Solution1

根据题目抓住关键点：

• 本题需要统计满足d(a) == b and d(b) == a, while a != b条件的数字对(a, b) ，故需要关注d(a) == a的情况；
• 统计因子之和的函数中，引入set数据结构，在结果中排除重复添加的情况，最后加上1。

def getSumOfDivisors(input_number):
    if input_number == 1 or input_number == 0:
        return 0
    important_set = set()
    for item in range(2, int(input_number/2)):
        if input_number % item == 0:
            important_set.add(item)
            important_set.add(input_number/item)
    return sum(important_set)+1

if __name__ == "__main__":
    blank_list = []
    for index_number in xrange(10000):
        if index_number in blank_list:
            continue
        tmp_result = getSumOfDivisors(index_number)
        if tmp_result == 0:
            continue
        if (index_number == getSumOfDivisors(tmp_result)) and (index_number != tmp_result) :
            blank_list.append(tmp_result)
            blank_list.append(index_number)
    print sum(blank_list)

Solution2

从第十题的分析可以看出，任何一个数字都可以表示成素数幂乘积的形式。下面对素数幂的因子之和进行推导分析：

设p为任意素数，由于素数只包含1和其本身两个因子，则：$\sigma(p) = p + 1$；

考虑p的a次幂，$\sigma(p^{a}) = 1 + p + p^{2} + p^{3} + … + p^{a}$ (1)；

式子(1)两侧同乘p，有$p*\sigma(p^{a}) = p + p^{2} + p^{3} + p^{4} + … + p^{a+1} (2)$;

(2) - (1)得，$p\sigma(p^{a}) - \sigma(p^{a}) = (p-1) \sigma(p^{a}) = p^{a+1} - 1$；

因此$\sigma(p^{a}) = (p^{a+1} - 1) / (p - 1)$

使用该公式并结合数字的素数幂分解即可较快求出因子之和，代码实现如下：

def SumOfDivisors(n):
    sum_result = 1
    p = 2
    # Prevents from checking prime factors greater than sqrt(n).
    while p**2 <= n and n > 1:
        if n % p == 0:
            j = p**2
            n /= p
            while n % p == 0:
                j *= p
                n /= p
            sum_result *= (j-1)
            sum_result /= (p-1)
        if p == 2:
            p = 3
        else:
            p += 2
    # Covers the case that one prime factor greater than sqrt(n) remains.
    if n > 1:
        sum_result *= (n+1)
    return sum_result

def SumOfProperDivisors(input_number):
    return SumOfDivisors(input_number) - input_number

与解法一对比，该方法将运算时间 从3s降低到0.3s。

22. Names score

Problem Description

Using names.txt (right click and ‘Save Link/Target As…’), a 46K text file containing over five-thousand first names, begin by sorting it into alphabetical order. Then working out the alphabetical value for each name, multiply this value by its alphabetical position in the list to obtain a name score.

For example, when the list is sorted into alphabetical order, COLIN, which is worth 3 + 15 + 12 + 9 + 14 = 53, is the 938th name in the list. So, COLIN would obtain a score of 938 × 53 = 49714.

What is the total of all the name scores in the file?

Solution

# Count the total scores of names
def SortContentAsAlphabeticalPosition(txtlocation):
    file = open(txtlocation, "r")
    content_string = file.read()
    content_list = content_string.replace("\"", "").split(",")
    return sorted(content_list)

def GetWordScore(input_word):
    word_Score = 0
    for item in input_word:
        word_Score += ord(item) - 65 + 1
    return word_Score

def CountTotalScore(content_list):
    total_score = 0
    for index in xrange(len(content_list)):
        total_score += GetWordScore(content_list[index]) * (index+1)
    return total_score
    
if __name__ == "__main__":
    txt_location = given_txt_path
    Sorted_list = SortContentAsAlphabeticalPosition(txt_location)
    total_score = CountTotalScore(Sorted_list)
    print total_score

23. Non-abundant sums

Problem Description

A perfect number is a number for which the sum of its proper divisors is exactly equal to the number. For example, the sum of the proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.

A number n is called deficient if the sum of its proper divisors is less than n and it is called abundant if this sum exceeds n.

As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number that can be written as the sum of two abundant numbers is 24. By mathematical analysis, it can be shown that all integers greater than 28123 can be written as the sum of two abundant numbers. However, this upper limit cannot be reduced any further by analysis even though it is known that the greatest number that cannot be expressed as the sum of two abundant numbers is less than this limit.

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

Solution

# Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers
def SumOfDivisors(n):
    sum_result = 1
    p = 2
    while p**2 <= n and n > 1:
        if n % p == 0:
            j = p**2
            n /= p
            while n % p == 0:
                j *= p
                n /= p
            sum_result *= (j-1)
            sum_result /= (p-1)
        if p == 2:
            p = 3
        else:
            p += 2
    if n > 1:
        sum_result *= (n+1)
    return sum_result
def SumOfProperDivisors(input_number):
    return SumOfDivisors(input_number) - input_number

if __name__ == "__main__":
    final_sum = 28123*28124/2
    least_abundant_number = 12
    possible_abundant_number_list = []
    while least_abundant_number <= 28123:
        for item in possible_abundant_number_list:
            if least_abundant_number - item in possible_abundant_number_list:
                final_sum -= least_abundant_number
                break
        if SumOfProperDivisors(least_abundant_number) > least_abundant_number:
            possible_abundant_number_list.append(least_abundant_number)
        least_abundant_number += 1
    print final_sum

使用上述暴力遍历的方法统计结果需要的时间较长，大概为160-180s左右。

24. Lexicographic permutations

Problem Description

A permutation is an ordered arrangement of objects. For example, 3124 is one possible permutation of the digits 1, 2, 3 and 4. If all of the permutations are listed numerically or alphabetically, we call it lexicographic order. The lexicographic permutations of 0, 1 and 2 are:

012 021 102 120 201 210

What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?

Solution

偷懒调用了python的itertools库，其中包含permutations()即排列函数。

# What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9
from itertools import permutations

if __name__ == "__main__":
    result_string = ""
    result_list = list(permutations('0123456789'))
    final_result = sorted(result_list)[999999]
    for item in final_result:
        result_string += item
    print eval(result_string)

25. 1000-digit Fibonacci number

Problem Description

The Fibonacci sequence is defined by the recurrence relation:

$$F_{n} = F_{n−1} + F_{n−2}, where F_{1} = 1 and F_{2} = 1$$

Hence the first 12 terms will be:

$$F_{1} = 1$$ $$F_{2} = 1$$ $$F_{3} = 2$$ $$F_{4} = 3$$ $$F_{5} = 5$$ $$F_{6} = 8$$ $$F_{7} = 13$$ $$F_{8} = 21$$ $$F_{9} = 34$$ $$F_{10} = 55$$ $$F_{11} = 89$$ $$F_{12} = 144$$

The 12th term, $F_{12}$, is the first term to contain three digits.

What is the index of the first term in the Fibonacci sequence to contain 1000 digits?

Solution

# What is the index of the first term in the Fibonacci sequence to contain 1000 digits

def GetIndexOfFibonacci():
    f = 1
    g = 1
    g_count = 2
    while len(str(g)) < 1000:
        f, g = g, f+g
        g_count += 1
    return g_count

if __name__ == "__main__":
    final_result = GetIndexOfFibonacci()
    print final_result

26. Reciprocal cycles

Problem Description

A unit fraction contains 1 in the numerator. The decimal representation of the unit fractions with denominators 2 to 10 are given:

1/2	=	0.5
1/3	=	0.(3)
1/4	=	0.25
1/5	=	0.2
1/6	=	0.1(6)
1/7	=	0.(142857)
1/8	=	0.125
1/9	=	0.(1)
1/10	=	0.1

Where 0.1(6) means 0.166666…, and has a 1-digit recurring cycle. It can be seen that 1/7 has a 6-digit recurring cycle.

Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.

Solution

本题为经典的循环节计算问题。

设分母为n，若n为合数，则必可以表示成多个素数的乘积，其循环节与这些分解素数中最大的循环节保持一致。故只需要讨论所给范围内所有素数分母其倒数的循环节即可。值得注意的是，1/2与1/5都是有限小数，故2和5这两个素数需要排除。

确定循环的范围后下面进行循环节计算：通过查阅资料可知，循环节问题可以等价为大整数分解问题。给定大整数n，求使得$10^{k} \equiv 1 (mod n)$成立的最小的k，该数值即为1/n循环节的长度。根据以上分析采取10**k-1 % n == 0作为限制条件求解。代码如下：

# Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part
# Brute force to find the longest recurring cycle
import math
# Judge whether n is prime
def isPrime(n):
    if n == 1:
        return False
    if n == 2:
        return True
    for i in range(2, int(math.sqrt(n))+1):
        if n % i == 0:
            return False
    return True

def GetRecurringCycle(input_number):
    L = 1
    while True:
        if (10**L-1) % input_number == 0:
            break
        L += 1
    return L

if __name__ == "__main__":
    max_recurring_length = 1
    result_n = 2
    for i in range(2, 1001):
        if isPrime(i) and i != 2 and i != 5:
            recurring_length = GetRecurringCycle(i)
            if recurring_length > max_recurring_length:
                max_recurring_length = recurring_length
                result_n = i
    print (result_n, max_recurring_length)

Mathematical Analysis

考虑上面提到的一般化问题：求最小的k使得$a^{k} \equiv 1(mod n)$，若n与a互素，求出分母n的欧拉函数值$\Phi(n)$，则循环节长度k必为其约数；若n与a存在公因子则无解。由此可见通过暴力遍历的方法可以求出约数k，从而得出循环节。

在RSA加密中，给定n=pq，则p，q均为大质数，此时1/n循环节的长度length为gcd(p-1, q-1)的约数。假定已知length的因数分解$length = l{1}^{c{1}}l{2}^{c{2}}l{3}^{c{3}}……*l{k}^{c{k}}$，则length共有$\prod[c_{i}+1]$个约数。将这些约数分别加上1，若某个约数y(j)加1后是质数，则y(j)+1可能是大整数n的约数。对所有小于$\sqrt{n}-1$的y(j)进行检验，必能找到一个恰好满足y(j)+1 = min(p, q)的数字。通过此种方法可以将大整数分解问题转换为求循环节的问题。在最坏的情况下，一个300位的大整数只需通过小于500次转换完成分解。

关于循环节长度计算公式存在如下改进：

\1.PNG)

该公式的应用如下：

\2.PNG)

参考论文：关于循环节长度计算公式的改进：https://www.ixueshu.com/document/6d924626ac37b4ee318947a18e7f9386.html

27. Quadratic primes

Problem Description

Euler discovered the remarkable quadratic formula:

$$n^{2} + n + 41$$

It turns out that the formula will produce 40 primes for the consecutive integer values 0≤n≤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.

The incredible formula $n^{2}−79*n+1601$ was discovered, which produces 80 primes for the consecutive values 0≤n≤79. The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form:

$$n^{2}+a*n+b , where|a|<1000 , |b|≤1000$$

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.

Solution

根据题目，可以仅通过数学推导的方法求解（条件极值问题+二次图像平移），参考链接如下：https://zhuanlan.zhihu.com/p/62137330。除此之外，还可以按照题意暴力求解，由于本题提供的a, b范围较小，故可以在较短时间内完成计算。根据题目给出的二次多项式$f(n) = n^{2} + a*n + b$可得出以下性质：

• 由f(0) = b可知参数b必为素数，故将b的遍历范围缩小为2—999；
• 由f(1) = a+b+1且该多项式为素数，可得结论a > -b-1，据此适当缩小a的遍历范围；
• 当gcd(a, b) = d且d != 1时，一定存在一个数n’ = minP(d)使得f(n’)为合数，其中minP(d)表示d的最小素因子。从而若a, b不互素，最多生成minP(d)个连续素数，从而在遍历过程中排除a, b不互素的情况即可。（需要用到math.gcd()函数且偷懒不想自己构造，故以下代码需选择Python3.5+版本运行）

遍历代码具有两个版本，记录如下：

Version1

# 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
# Some properties:
    # f(0) = b, so b must be a prime
    # a > -b-1
    # if gcd(a, b) = d and d != 1, we could only get as many as the least prime factor of d continuous primes
import math

def isPrime(n):
    if n <= 0 or n == 1:
        return False
    if n == 2:
        return True
    for i in range(2, int(math.sqrt(n))+1):
        if n % i == 0:
            return False
    return True

def countPrimes(a, b):
    prime_counter = 0
    n = 1
    while True:
        if isPrime(n*n+a*n+b):
            prime_counter += 1
        else:
            break
        n += 1
    return prime_counter


if __name__ == '__main__':
    result_set = []
    max_result = 0
    for a in range(-995, 999, 2):
        for b in range(2, 1000):
            if isPrime(b) and math.gcd(a, b) == 1:
                count = countPrimes(a, b)
                if count > max_result:
                    max_result = count
                    result_set.append(a*b)
            else:
                continue
    print(result_set[-1])

Version2

引入python中的数学运算库sympy，该库继承多种数学运算，包括大量素数操作的接口：

isprime(n)              # Test if n is a prime number (True) or not (False).
 
primerange(a, b)        # Generate a list of all prime numbers in the range [a, b).
randprime(a, b)         # Return a random prime number in the range [a, b).
primepi(n)              # Return the number of prime numbers less than or equal to n.
 
prime(nth)              # Return the nth prime, with the primes indexed as prime(1) = 2. The nth prime is approximately n*log(n) and can never be larger than 2**n.
prevprime(n, ith=1)     # Return the largest prime smaller than n
nextprime(n)            # Return the ith prime greater than n
 
sieve.primerange(a, b)  # Generate all prime numbers in the range [a, b), implemented as a dynamically growing sieve of Eratosthenes.

代码如下，其中应用到sympy.isprime()函数进行素数判断，应用sieve.primerange(a, b)函数生成1—1000内的所有素数。

import sympy

# This function is used to collect all the primes between 1 and 1000
Prime = list(sympy.sieve.primerange(1, 1000))
dic = {}
dic['n'] = 0
for a in range(-995, 999, 2):
    for b in Prime:
        n = 0
        c = n*n + a*n + b
        while sympy.isprime(c):
            n += 1
            c = n*n + a*n + b
        n -= 1
        if n > dic['n']:
            dic['n'] = n
            dic['a*b'] = a*b
print dic['a*b']

sympy库函数文档：http://docs.sympy.org/latest/modules/ntheory.html?highlight=prime#ntheory-functions-reference

28. Number spiral diagonals

Problem Description

Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:

\3.PNG)

It can be verified that the sum of the numbers on the diagonals is 101.

What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral formed in the same way?

Solution

# gridA: 5*5 there are 3 circles
# gridB: 7*7 there are 4 circles
# gridN: N*N there must be (N+1)/2 circles
# Arithmetic progression

def getCircles(n):
    return (n+1) / 2

def countSum(circle_number):
    final_sum = 0
    step = 2
    origin_number = 3
    for i in range(circle_number-1):
        final_sum += 4*origin_number + 6*step
        origin_number += 4*step + 2
        step += 2
    return final_sum
if __name__ == '__main__':
    final_result = countSum(getCircles(1001))
    print final_result+1

29. Distinct powers

Problem Description

Consider all integer combinations of $a^{b}$ for 2 ≤ a ≤ 5 and 2 ≤ b ≤ 5:

$$2^{2}=4, 2^{3}=8, 2^{4}=16, 2^{5}=32$$ $$3^{2}=9, 3^{3}=27, 3^{4}=81, 3^{5}=243$$ $$4^{2}=16, 4^{3}=64, 4^{4}=256, 4^{5}=1024$$ $$5^{2}=25, 5^{3}=125, 5^{4}=625, 5^{5}=3125$$

If they are then placed in numerical order, with any repeats removed, we get the following sequence of 15 distinct terms:

4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125

How many distinct terms are in the sequence generated by $a^{b}$ for 2 ≤ a ≤ 100 and 2 ≤ b ≤ 100?

Solution

final_set = set()
for a in range(2, 101):
	for b in range(2, 101):
		tmp_a_b = a**b
        tmp_b_a = b**a
        if tmp_a_b not in final_set:
            final_set.add(tmp_a_b)
        if tmp_b_a not in final_set:
            final_set.add(tmp_b_a)
print len(final_set)

30. Digit fifth powers

Problem Description

Surprisingly there are only three numbers that can be written as the sum of fourth powers of their digits:

$$1634 = 1^{4} + 6^{4} + 3^{4} + 4^{4}$$ $$8208 = 8^{4} + 2^{4} + 0^{4} + 8^{4}$$ $$9474 = 9^{4} + 4^{4} + 7^{4} + 4^{4}$$

As$ 1 = 1^{4} $is not a sum it is not included.

The sum of these numbers is 1634 + 8208 + 9474 = 19316.

Find the sum of all the numbers that can be written as the sum of fifth powers of their digits.

Solution

该问题的关键为确定枚举的上界：设满足题设条件的数字为n位数，取其每位都为9的情况，构造不等式$ 10^{n} <= 9^{5}*n $，计算满足条件n的最小值即可找出上界。代码如下：

# Find the sum of all the numbers that can be written as the sum of fifth powers of their powers
def getSum(input_number):
    sum_result = 0
    while input_number != 0:
        sum_result += (input_number % 10)**5
        input_number /= 10
    return sum_result

if __name__ == '__main__':
    final_result = 0
    # Work out the upper bound of the enumeration range
    # 10^n <= 9^5*n => we can get the minimum value of n is 5
    for item in range(2, 9**5*6):
        if item == getSum(item):
            final_result += item
            print item
    print final_result