欧拉计划25-40

10 minute read

Published: March 16, 2022

所需包

import math
import time
import numpy as np

第二十五题

这是一个数列递推问题，如果使用递归会使用O($n^2$)的空间与O($n*2^n$)的时间，我们直接使用递推公式一步一步推

start = time.time()

n = 1000
m = 10**(n - 1)
a_1 = 1
a_2 = 1
ans = 2
while a_2 < m:
    ans += 1
    a_1, a_2 = a_2, a_1 + a_2
print("第一个超过{}位的数字是位于第{}位的:\n{}".format(n, ans, a_2))


print("消耗时间{:.6f}s".format(time.time() - start))

第一个超过1000位的数字是位于第4782位的:
1070066266382758936764980584457396885083683896632151665013235203375314520604694040621889147582489792657804694888177591957484336466672569959512996030461262748092482186144069433051234774442750273781753087579391666192149259186759553966422837148943113074699503439547001985432609723067290192870526447243726117715821825548491120525013201478612965931381792235559657452039506137551467837543229119602129934048260706175397706847068202895486902666185435124521900369480641357447470911707619766945691070098024393439617474103736912503231365532164773697023167755051595173518460579954919410967778373229665796581646513903488154256310184224190259846088000110186255550245493937113651657039447629584714548523425950428582425306083544435428212611008992863795048006894330309773217834864543113205765659868456288616808718693835297350643986297640660000723562917905207051164077614812491885830945940566688339109350944456576357666151619317753792891661581327159616877487983821820492520348473874384736771934512787029218636250627816
消耗时间0.001000s

第二十六题

经过手动列式除法可见，当余数是之前出现过的时候，陷入循环

def reciprocal_cycles(x: int)->int:
    if x <= 0: return -1
    digits = []#记录曾经出现过的余数
    n = 1
    a = n%x
    while a not in digits:
        digits.append(a)
        n *= 10
        a = n%x
    if a == 0: return 0
    return (len(digits) - digits.index(a))
    
start = time.time()

n = 1000
ans = 0
v = -1
for d in range(n):
    a = reciprocal_cycles(d)
    if a > v:
        ans = d
        v = a
print("{}拥有最长的循环列，长度为{}".format(ans, v))

print("消耗时间{:.6f}s".format(time.time() - start))

983拥有最长的循环列，长度为982
消耗时间0.212047s

第二十七题

很明显我们需要多次使用验证质数，为了节省时间，先用欧拉筛找出素数

def prime(n: int) -> list:
    '''
    :param n: 限制n
    :return: 不超过n的质数全体
    '''
    n = n + 1
    num = [True] * n
    num[0:2] = [False] * 2
    for i in range(n):
        if num[i]:
            j = 2
            while i * j < n:
                num[i * j] = False
                j += 1
    return [i for i in range(0, n) if num[i]]

start = time.time()

primes = prime(1000000)

print("消耗时间{:.6f}s".format(time.time() - start))

消耗时间0.534806s

$n=0$时，多项式取$b$的值

这意味着，最长序列的b必须是质数，否则用题干给的多项式会产生更长的序列

接着挨个验证, 可能需要花上一点时间（大约5分钟）

def quadratic_primes(a, b):
    ans = 0
    while ans**2 + a* ans + b in primes:
        ans += 1
    return ans


start = time.time()

n_a = 1000
n_b = 1000

max_l = 0
ans = 0
for b in [p for p in primes if p <= n_b]: # 当n = 0 时，多项式取值为b
    for a in range(-n_a-1, n_a):
        if 1 + a + b in primes: # 再减少一部分需要验证的(好像没有优化多少)
            c = quadratic_primes(a,b)
            if c>max_l:
                max_l = c
                ans = a*b
print(ans, max_l)

print("消耗时间{:.6f}s".format(time.time() - start))

-59231 71
消耗时间328.339104s

第二十八题

法一，暴力至上先得到相应的矩阵，再累加

空间：$O(n^2)$ , 时间：$O(n^2)$

start = time.time()

n = 1001
grid = [[0]*n for i in range(n)]

# 制作符合要求的矩阵
i, j = n // 2, n // 2
grid[i][j] = 1
c = 1
for a in range(n//2):
    c += 1
    j += 1
    grid[i][j] = c
    for b in range(2*(a+1)-1):
        i += 1
        c+=1
        grid[i][j] = c
    for b in range(2*(a+1)):
        j -= 1
        c+=1
        grid[i][j] = c
    for b in range(2*(a+1)):
        i -= 1
        c+=1
        grid[i][j] = c
    for b in range(2*(a+1)):
        j += 1
        c+=1
        grid[i][j] = c

# 累加
ans = 0
for i in range(n):
    ans += grid[i][i]
    ans += grid[i][n-1-i]
ans -= 1
print("{}*{}螺旋数阵对角线元素之和为{}".format(n, n, ans))

print("消耗时间{:.6f}s".format(time.time() - start))

1001*1001螺旋数阵对角线元素之和为669171001
消耗时间0.158036s

或者，简单分析一下，第一圈有1个数，第二圈有8个数，第三圈有16个数, …, 第n圈有8*(n-1)个

是$(2n-1)*(2n-1)$矩阵

第$n$圈增加的数值恰好是$((2n-3)^2 + 12(n-1))+((2n-3)^2 + 22(n-1)))+((2n-3)^2 + 23(n-1))+((2n-3)^2 + 24(n-1))$

整理一下每圈增加$4(2n-3)^2 + 20(n-1) = 16n^2 - 28n + 16$

再进行累加 $\begin{array}{l} 1 + \sum_{k = 2}^{n} (16k^2 - 28k + 16) \\ = 1 + 16*\frac{n*(2n+1)*(n+1)}{6} - 14*n(n+1) + 16*n - 4\\ = \frac{16}{3} n^3 - 6 n^2 + \frac{14}{3} n - 3 \end{array}$

这下时间空间都是$O(1)$

start = time.time()

n=1001
k = n//2 + 1
print((16/3)*k**3 - 6*k**2 + (14/3)*k -3) 

print("消耗时间{:.10f}s".format(time.time() - start))

669171001.0
消耗时间0.0000000000s

第二十九题

暴力计数, 几乎花不了多少时间

start = time.time()

n_a = range(2,101)
n_b = range(2,101)

nums = []
for a in n_a:
    nums += [a**b for b in n_b]
print(len(set(nums)))

print("消耗时间{:.6f}s".format(time.time() - start))

9183
消耗时间0.005002s

第三十题

先找出这个问题的上限

n = 5
ans = 1
while (10**ans - 1) < ans * (9**n):
    ans += 1
print("最大能满足条件的数字不超过{}".format(ans*(9**n)))

最大能满足条件的数字不超过354294

然后暴力（反正花不了多少时间

start = time.time()

n = 5
num = []
for a in range(0, 4):
    for b in range(0, 10):
        for c in range(0,10):
            for d in range(0,10):            
                for e in range(0,10):
                    for f in range(0,10):
                        m = a*10**5 + b*10**4 + c*10**3 + d*10**2 + e*10 + f
                        if m > 354294: continue # 跳过明显不满足条件的数
                        p = a**n + b**n + c**n + d**n + e**n + f**n
                        if m == p:
                            num.append(m)
print(num)
print(sum(num[2:])) # 消除0，1这两个平凡解

print("消耗时间{:.6f}s".format(time.time() - start))

[0, 1, 4150, 4151, 54748, 92727, 93084, 194979]
443839
消耗时间0.843487s

第三十一题

暴力枚举，当计算传递到1时，无论还剩多少都只有一种补全方式,计数，

PS：看起来一点都不具有美感

start = time.time()

POUND2=200
def countP100(need):
    p = 100
    count = 0
    for c in range(0,need//p+1):
        count += countP50(need - p*c)
    return count
def countP50(need):
    p = 50
    count = 0
    for c in range(0,need//p+1):
        count += countP20(need - p*c)
    return count
def countP20(need):
    p = 20
    count = 0
    for c in range(0,need//p+1):
        count += countP10(need - p*c)
    return count
def countP10(need):
    p = 10
    count = 0
    for c in range(0,need//p+1):
        count += countP5(need - p*c)
    return count
def countP5(need):
    p = 5
    count = 0
    for c in range(0,need//p+1):
        count += countP2(need - p*c)
    return count
def countP2(need):
    p = 2
    count = 0
    for c in range(0,need//p+1):
        count += 1
    return count


print(countP100(POUND2)+1) # 记得计数P200 本身

print("消耗时间{:.6f}s".format(time.time() - start))

73682
消耗时间0.003001s

稍微整理一下思路，或许我们可以用以下等式表示我们的操作 $w(t,c) = \begin{cases} 1 & \text{ if } c=1 \vee t=0 \\ \sum_{i=0}^{\lfloor \frac{t}{c} \rfloor } w(t-i*c, s(c)) & \text{ if } c>1 \wedge t>0 \end{cases}$ 其中，$t$代表需要的金额，$c$代表使用的钱币，$s(c)$将给出面值序列中$c$的下一个值

借着这个递归式，我们可以写出泛化能力更强的算法（PS：也更优雅

start = time.time()

POUND2=200
cents = [200,100,50,20,10,5,2,1]

def w(t: int, n: int=0)->int:
    c = cents[n]
    if c == 1 or t == 0: return 1
    count = 0
    for i in range(0, t//c + 1): 
        count += w(t-i*c, n+1)
    return count

print(w(POUND2))

print("消耗时间{:.6f}s".format(time.time() - start))

73682
消耗时间0.021535s

第三十二题

每日暴力验证（1/1

至于上限为什么是$10^4$,那是因为最差的部分大概就是 4位数×1位数=4位数

start = time.time()

ans = set()
for a in range(10000):
    for b in range(a):
        num_t = str(a) + str(b) + str(a*b)
        if len(num_t) != 9: continue
        for c in range(1,10):
            if str(c) not in num_t: break
            if c == 9 and str(c) in num_t: ans.add(a*b)
        
print(len(set(ans)))
print(sum(set(ans)))

print("消耗时间{:.6f}s".format(time.time() - start))

7
45228
消耗时间40.450104s

但有些东西似乎可以优化，比如对于一个a而言，num_t长度过长的时候，其余更大的b值都不需要验证

可以看到，效果拔群。几十秒经过优化变成几十毫秒

start = time.time()

ans = set()
for a in range(10000):
    for b in range(a):
        num_t = str(a) + str(b) + str(a*b)
        if len(num_t) < 9: continue
        if len(num_t) > 9: break
        for c in range(1,10):
            if str(c) not in num_t: break
            if c == 9 and '9' in num_t: ans.add(a*b)
        
# print(len(set(ans)))
print(sum(ans))

print("消耗时间{:.6f}s".format(time.time() - start))

45228
消耗时间0.099533s

第三十三题

先用辗转相减法定义最大公因数函数

def gcf(a: int, b: int)->int:
    if a == b: return a
    if a < b: a, b = b, a
    return gcf(a - b, b)
gcf(132,111)

对题干进行分析可知相当于在解这样一个方程 $a \% 10 == b // 10$ 分数值小于1相当于 $a < b$ 这样可以轻松搜出所有解

start = time.time()

p = 1
q = 1
for a in range(11, 100):
    for b in range(a+1, 100):
        if b//10 == a%10 and b%10 * a == a//10 * b:
            p *= a
            q *= b

c = gcf(p,q)
print(q / c)


print("消耗时间{:.6f}s".format(time.time() - start))

100.0
消耗时间0.002001s

第三十四题

def factorial(n: int)->int:
    ans = 1
    for i in range(2, n+1):
        ans *= i
    return ans

遇到没给范围的题，先搜上限

显然，我们有

ans = 1
while 10**ans - 1 < ans* factorial(9):
    ans += 1
print("解空间将不会超过{}".format(ans* factorial(9)))

解空间将不会超过2540160

start = time.time()

factorials = [factorial(x) for x in range(10)] # 提前算好，防止无效重复
func = lambda x: sum([factorials[int(t)] for t in str(x)])

ans = 0
for i in range(3, 2540160):
    if func(i) == i:
        ans+=i

print(ans)
        
print("消耗时间{:.6f}s".format(time.time() - start))

40730
消耗时间6.404168s

第三十五题

先引进欧拉筛，因为需要大量判别，所以我们这次引用一个判别数组

def prime(n: int) -> list:
    '''
    :param n: 限制n
    :return: 一个代表是否是质数的bool数组
    '''
    n = n + 1
    num = [True] * n
    num[0:2] = [False] * 2
    for i in range(n):
        if num[i]:
            j = 2
            while i * j < n:
                num[i * j] = False
                j += 1
    return num

start = time.time()

primes = prime(1000000)

print("消耗时间{:.6f}s".format(time.time() - start))

消耗时间0.507079s

循环检测，如果有一个不是质数返回False，如果全是质数返回True 然后找到所有符合条件的数然后计数就好

def is_circle_prime(x: int):
    num_t = str(x)
    for i in range(0, len(num_t)):
        if not primes[int(num_t)]:
            return False
        num_t = num_t[1:] + num_t[0]
    return True

start = time.time()

primes = prime(1000000)
print(len([i for i in range(1000000) if is_circle_prime(i)]))

print("消耗时间{:.6f}s".format(time.time() - start))

55
消耗时间0.906713s

第三十六题

暴力验证

因为不考虑前导零，所以该数字必然是奇数（否则二进制最后一位是0

利用内置函数将n分别变成10进制表示和2进制表示，然后验证是否可反转

def is_palindrome(t: str)->bool:
    return t == t[::-1]

start = time.time()

ans = sum([i for i in range(1, 1000000, 2) if is_palindrome(str(i)) and is_palindrome(str(bin(i))[2:])])
print(ans)

print("消耗时间{:.6f}s".format(time.time() - start))

872187
消耗时间0.165039s

欧拉计划法一

尝试欧拉计划给出的文档中的算法

考虑n进制，一位一位的将反向序列算出来，相比于我的方法而言，它的独立性更高

def isPalindrome(x: int, base: int)->bool:
    reversed = 0
    k = x
    while k > 0:
        reversed = reversed*base + k%base
        k //= base
    return reversed == x

start = time.time()

ans = sum([i for i in range(1, 1000000, 2) if isPalindrome(i, 10) and isPalindrome(i, 2)])
print(ans)

print("消耗时间{:.6f}s".format(time.time() - start))

872187
消耗时间0.613241s

第三十七题 Truncatable_primes

一个一个测试，如果有某个数量级不存在可截素数，那么，以后也不会有，借此我们获得了停止条件

def prime(n: int) -> list:
    '''
    :param n: 限制n
    :return: 一个代表是否是质数的bool数组
    '''
    n = n + 1
    num = [True] * n
    num[0:2] = [False] * 2
    for i in range(n):
        if num[i]:
            j = 2
            while i * j < n:
                num[i * j] = False
                j += 1
    return num

start = time.time()

limit = 10**7
primes = prime(limit)

print("消耗时间{:.6f}s".format(time.time() - start))

消耗时间3.378602s

def is_left_truncatable_ptime(x: int)->bool:
    if not primes[x]: return False
    if x < 10: return True
    return is_left_truncatable_ptime(int(str(x)[1:]))

def is_right_truncatable_prime(x: int)->bool:
    if not primes[x]: return False
    if x < 10: return True
    return is_right_truncatable_prime(x//10)


start = time.time()

ans = []
for i in range(1, 7):
    ok = False
    for j in range(10**i, 10**(i+1)):
        if not ok and (is_right_truncatable_prime(j) or is_left_truncatable_ptime(j)): ok = True
        if is_left_truncatable_ptime(j) and is_right_truncatable_prime(j):
            ans.append(j)
    if not ok: break

print(ans)
print(sum(ans))

print("消耗时间{:.6f}s".format(time.time() - start))

[23, 37, 53, 73, 313, 317, 373, 797, 3137, 3797, 739397]
748317
消耗时间1.322324s

需要花上一点时间，可以想想有没有什么优化的方法

第三十八题Pandigital multiples

考虑一个数x，如果它有5位，则与1，2分别相乘后则至少有10位数字

显然，$10^4$是这个问题的一个上届

def pandigital_multiples(x: int):
    nums = [False]*10
    ans = ''
    i = 1
    while True:
        a = x*i
        while a > 0:
            b = a%10
            if b == 0 or nums[b]: return (False, 0)
            nums[b] = True
            a //= 10
        if all(nums[1:]): return (True, int(''.join([str(x*n) for n in range(1, i+1)])))
        i += 1

start = time.time()

ans = 0
for i in range(1, 10**5):
    (ok, num) = pandigital_multiples(i)
    if ok and num > ans:
        ans = num
print(ans)
        
print("消耗时间{:.6f}s".format(time.time() - start))

932718654
消耗时间0.099027s

第三十九题

取出每一个进行验证即可

start = time.time()

p = 1000
nums = [0]*(p+1)
for a in range(1, p):
    for b in range(a, p):
        c = math.sqrt(a**2 + b**2)
        if int(c) == c and a + b + c <= p:
            nums[a+b+int(c)] += 1

ans = nums.index(max(nums))
print(ans)

print("消耗时间{:.6f}s".format(time.time() - start))

840
消耗时间0.425792s

第四十题 Champernowne’s constant

因为题干中出现了$d_{1000000}$，如果直接将数字打印出来，则会产生大约1M位的数字。

显然，我们其实并不需要这么多的空间，因为计算$d_{1000000}$的时候我们并不需要知道$d_1 - d_{999999}$的值，所以，我们考虑通过索引直接计算数字的办法

先找一下规律，一位数有9个，两位数有90个，三位数有900个，。。。，n位数有$9*10^{n-1}$个

对应到数字数，一位数有$91$个，两位数有$902$个，三位数有$9003$个，。。。，n位数有$910^{n-1}*n$个

那么我们可以用一个循环确定位数，当确定位数之后，再定位到第几个数的第几位。

def champernowne_constant(index: int)->int:
    digits = 1
    while index >= digits*9*10**(digits-1):
        index -= digits*9*10**(digits-1)
        digits += 1

    num = 10**(digits-1) + index//digits

    for i in range(digits - index%digits - 1):
        num //= 10
    return num%10


start = time.time()

# d1 × d10 × d100 × d1000 × d10000 × d100000 × d1000000
nums = [1,10,100,1000,10000,100000,1000000]
ans = 1
for d in nums:
    ans *= champernowne_constant(d-1)# 这里使用 d-1 是因为我的索引是从 0 开始的
print(ans)

print("消耗时间{:.6f}s".format(time.time() - start))

210
消耗时间0.000000s

Share on

Twitter Facebook LinkedIn

Zijie Wen