欧拉计划41-50
Published:
import time
import math
第四十一问 全数字的素数
这个问题有一个显然上限 987654321, 但是对于我们的小计算机来说，它太大了。但是，我们先验证是否是全序列花费的时间远比验证是否是质数的时间短，这样我们可以用尽可能少的算力排除大量的不可能事件
因为只需要验证是否是质数，那么我们其实只要验证是否能被$1 - \sqrt{987654321}$的质数整除即可
首先邀请我们的老朋友——欧拉筛，入场，计算出$1 - \sqrt{987654321}$的质数序列用来验证
def prime(n: int) -> list:
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(int(math.sqrt(987654321)))
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
你浪费了12.00ms
然后，我们定义一个验证是否是质数的函数
def is_prime(x):
limit = math.sqrt(x)
for p in primes:
if x % p == 0:
return False
if p > limit:
break
return True
至于是不是全序列，不需要测量，我们直接用排列组合枚举全序列数，这样我们需要测试的数字只有$9! = 362880$个远远少于$987654321$
from itertools import permutations # 懒得重复造轮子了
def pandigital_prime():
nums = '987654321'
for i in range(9):
for x in permutations(nums[i:]): # 因为函数特性，将从最大数开始枚举
n = int(''.join(x))
if is_prime(n):
return n
start = time.time()
print(pandigital_prime())
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
7652413
你浪费了238.13ms
第四十二问 编码三角形数
题目给的数据是一个txt文件，我们需要先将他读取进来并处理成列表形式
with open('source/p042_words.txt', encoding='utf-8') as f:
words = f.read().replace('\"','').split(',')
f.close()
max_l = max([len(w) for w in words])
print(max_l)
14
我们注意到，这个单词列表里最长的单词也只有14个单词的长度，也就是说，单词编码不会超过$14 \times 26 = 364$
单词不多，只有两千个，我们直接遍历就好
start = time.time()
def word_value(word):
return sum([ord(t) - 64 for t in word])
tri_nums = [n*(n+1)//2 for n in range(1, int(math.sqrt(max_l*26*2)) +1)]
print(len( [w for w in words if word_value(w) in tri_nums]))
print("你浪费了{:.6f}ms".format((time.time()-start)*1000))
162
你浪费了1.000166ms
第四十三问 子串的可整除性
from itertools import permutations
start = time.time()
divisors = [2,3,5,7,11,13,17]
ans = 0
for num in permutations(list(range(0,10))):
if num[0] == 0: continue
ok = True
for i in range(7):
n = num[i+1]*100 + num[i+2]*10 + num[i+3]
if n % divisors[i] != 0:
ok = False
break
if ok:
ans += int(''.join([str(k) for k in num]))
print(ans)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
16695334890
你浪费了2250.26ms
第四十四问 五边形数
暂时没想到好办法，所以直接暴力了
什么？你问我上界怎么来的？试出来的呗.
start = time.time()
p=[i*(3*i-1)//2 for i in range(1,5000)]
k=[False]*(5000*(3*5000-1))
for i in p:
k[i]=True
ans=190000000
for i in p:
for j in p:
if k[abs(i-j)] and k[i+j]:
# print(i, j)
ans = min(ans,abs(i-j))
print(ans)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
5482660
你浪费了6323.92ms
第四十五问 三角形数、五边形数和六角形数
这一次，我们试一下用生成器
一个一个数生成，不断更新三个数列中的最小值，直到三个数列取到相同的数字
from itertools import count
def triangular():
for n in count(1):
yield n*(n+1)//2
def pentagonal():
for n in count(1):
yield n*(3*n-1)//2
def hexagonal():
for n in count(1):
yield n*(2*n-1)
def nums():
a = triangular()
b = pentagonal()
c = hexagonal()
while True:
x = next(a)
y = next(b)
z = next(c)
while True:
if x == y and y == z:
break
if x > y:
y = next(b)
if y > x:
x = next(a)
if y > z:
z = next(c)
if y < z:
y = next(b)
yield x
start = time.time()
a = nums()
for i in range(3):
print(next(a))
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
1
40755
1533776805
你浪费了32.01ms
更进一步，我们发现第一个序列包含了第三个序列，这意味着我们只需要验证后两个序列即可
from itertools import count
def hexagonal():
for n in count(1):
yield n*(2*n-1)
def pentagonal():
for n in count(1):
yield n*(3*n-1)//2
def nums():
a = hexagonal()
b = pentagonal()
while True:
x = next(a)
y = next(b)
while True:
if x == y:
break
if x > y:
y = next(b)
if y > x:
x = next(a)
yield x
start = time.time()
a = nums()
for i in range(3):
print(next(a))
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
1
40755
1533776805
你浪费了15.99ms
第四十六问 哥德巴赫的另一个猜想
或许我们可以参照欧拉筛的思想做一个 “哥德巴赫的另一个猜想”筛（非专业术语，我瞎叫的
def prime(n: int) -> list:
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()
n = 10000
is_prime = prime(n)
primes = [i for i in range(0, n) if is_prime[i]]
squares = [i**2 for i in range(1, int(math.sqrt(n)))]
nums = [False]*(n+1)
for p in primes:
for q in squares:
if p+2*q > n: break
nums[p+2*q] = True
ans = [i for i,y in enumerate(nums) if not y and not is_prime[i] and i&1 == 1]
print(ans)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
[1, 5777, 5993]
你浪费了14.00ms
第四十七问 不同的质因数
本题最大的关键在于如何确定一个数字的质因数的数量。反正不用管质因数的幂数，直接见一个就直接除尽一个。
以及为了迭代方便，我们使用双指针，因为当我们检测到一个数不符合要求的时候，那么任何一列包含这个数的数都不符合要求，我们直接跳到这个数后面继续就好，这样就不用对同一个数反复检测。
start = time.time()
def prime(n: int) -> list:
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]]
primes = prime(10**5)
def count_factor(x):
ans = 0
for p in primes:
if p > x: break
if x % p == 0:
ans += 1
while x%p == 0:
x = x//p
return ans
n = 4
i = 0
j = 0
while True:
if j - i == n:
break
if count_factor(j) == n:
j += 1
else:
j += 1
i = j
print([x for x in range(i,j)])
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
[134043, 134044, 134045, 134046]
你浪费了6954.51ms
第四十八问 自幂
取一个数字的后十位，相当于是对一个数取模，虽然说python不需要向c一样考虑整数溢出的事情，但为了效率，我们最好尽量避免处理太大的数，那样子会侵占过多的计算机资源，同时也会消耗过多的时间。
对于取模后的运算，我们可以验证以下定理成立
\(a \pmod{m} + b \pmod{m} \equiv (a+b) \pmod{m}\) \(a \pmod{m} * b \pmod{m} \equiv (a*b) \pmod{m}\)
这意味着，我们可以通过在计算加法与乘法之前取模来减少对不必要的数的处理
start = time.time()
n = 1000
mod = 10**10
def mod_pow(x, p, m):
ans = 1
for i in range(p):
ans *= x
ans %= m
return ans
ans = 0
for i in range(1, n+1):
ans += mod_pow(i, i, mod)
ans %= mod
print("{:0>10d}".format(ans))
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
9110846700
你浪费了55.57ms
第四十九问 素数重排
这一道题第一眼看上去的时候感觉有点无从下手
先试试暴力枚举
反正，运算的任务是芯片，储存的任务是内存，跟我又有什么关系呢
start = time.time()
length = 4
down = 10**(length-1)
def prime(n: int) -> list:
global down
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(down, n) if num[i]]
primes = prime(10**length)
for a in primes:
for b in primes:
if b >= a: break
for c in primes:
if c >= b: break
if (2*b == a+c) and (set(str(a)) == set(str(b)) == set(str(c))):
print(c,b,a)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
1487 4817 8147
2969 6299 9629
你浪费了21861.58ms
嘶，用的时间也太久了吧。
回头有时间优化一下
第五十问 连续素数的和
def prime(n: int) -> list:
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]], num
暴力验证（1/1 空间：$O(n)$ ，时间：$O(n^3)$ 消耗时间直接破出天际
start = time.time()
limited = 1000000
primes, _ = prime(limited)
n = len(primes)
def consecutive_prime_sum(x):
global primes
ans = -1
for i in range(n):
if primes[i] >= x: break
c = 0
v = 0
for j in range(i,n):
if v >= x: break
v += primes[j]
c += 1
if v == x:
ans = max(ans, c)
return ans
ans = 0
v = -1
for i in primes:
if consecutive_prime_sum(i) > v:
v = consecutive_prime_sum(i)
ans = i
print(ans, v)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
上一百万花的时间太久，感兴趣的可以用自己的电脑逝一逝
懒得算了，换个思路
因为是连续序列，我们可以先遍历素数求和，这样时间复杂度会缩到O(n^2),接着，其实不用每次都计算和，我们注意到每次只会增加一个素数，这样时间将进一步优化。
start = time.time()
limited = 1000000
primes, is_prime = prime(limited)
n = len(primes)
a,b = 0,0
ans = 0
v = -1
for a in range(n):
d = 0
for b in range(a,n):
d += primes[b]
if d > limited: break
if is_prime[d]:
if b-a > v:
v = b-a
ans = d
print(ans,v+1)
print("你浪费了{:.2f}ms".format((time.time()-start)*1000))
997651 543
你浪费了437.10ms
完成~(^o^)~~