欧拉计划

https://pe-cn.github.io/

”欧拉计划是一系列有挑战性的数学与计算机编程题；要解开它们，需要的不止是数学知识：尽管数学能够帮助你找到一些优雅而有效的方法，大多数题目仍需要借助计算机和编程技巧来完成解答。

创立欧拉计划的初衷，以及不断维持其运行的动力，在于为好奇的头脑提供一个平台，使他们能够在有趣愉悦的氛围中，探索未知领域，学习新的知识。

欧拉计划预期的受众，包括在基础课程外学有余力的学生、非数学背景但对数学感兴趣的成年人以及希望磨炼解题能力或是数学能力的专业人士。

欧拉计划的题目难度不一，对于大多数人来说，解题就是一个逐渐学习的过程。也就是说，每当你解开一个题目，你将会了解一些新的知识，从而帮助你着手解决以前无从下手的题目。因此，任何有决心的参与者，即使进展再缓慢，也一定能逐一解开每一道题。

欧拉计划的存在，是为了每个对数学感兴趣的人，鼓励他们，挑战他们，并最终培养他们的能力与乐趣。“

Problem 46 题目发布于 2003-06-21 翻译更新于 2015-07-31

克里斯蒂安·哥德巴赫曾经猜想，每个奇合数可以写成一个素数和一个平方的两倍之和。

• 9 = 7 + 2×1²
• 15 = 7 + 2×2²
• 21 = 3 + 2×3²
• 25 = 7 + 2×3²
• 27 = 19 + 2×2²
• 33 = 31 + 2×1²

最终这个猜想被推翻了。

最小的不能写成一个素数和一个平方的两倍之和的奇合数是多少？

答案可能是5777：

import math
primes = set()
N = 6000
for i in range(1, N):
	is_prime = True
	for j in range(2, i):
		if i % j == 0:
			is_prime = False
			break
	if is_prime: primes.add(i)
#print('"primes":', primes)
pl = list(sorted(primes))
#print(pl)

#for i in range(33, 5772, 2):
for i in range(5771, N, 2):
	if i in primes: continue
	good = True
	for j in pl:
		if i < j: break
		t = i - j
		if t % 2 == 0:
			t = math.sqrt(t / 2)
			#if t % 1 == 0:
			if int(t) == t:
				print(i, '=', j, '+2*', int(t), '^2', sep='')
				good = False; break
	if good: print('*****', i); break

for i in range(1, 500):
	for j in range(1, i):
		if i * j == 5777: print(i, j); quit()

如何判断一个数是不是平方数？好像用集合最快，我跑偏了一小会。搞个表在sqrt之前没啥意思。Icarus的sample里有算sqrt的: https://files.cnblogs.com/files/blogs/714801/sqrt.zip

python版: 面试题, floor(sqrt(x))的几种方法 - 博客园

# 土了土了，忘了著名的Quake公式(n+1)^2 - n^2 = 2n+1了，即循环步长为2n+1. 可设置100, 1万等几个检查站，
# 不必每次从1开始。不过2^17=65536*2，集合17步搜到。安排add的顺序使二叉树平衡。集合费内存。
# x^0.5一求导，变得好像越来越难算了。
import math

s = set()
for i in range(1<<17): s.add(i*i)
print(len(s))

from numba import jit
@jit(nopython=True)
def is_square_1(x):
	x = math.sqrt(x)
	return int(x) == x

#@jit(nopython=True)
def is_square_2(x): return x in s

@jit(nopython=True)
def is_square_3(x):
	i = 1
	while True:
		t = i * i
		if t == x: return True
		elif t > x: break
		i += 100
	i -= 100
	while True:
		t = i * i
		if t == x: return True
		elif t > x: break
		i += 10
	i -= 10
	while True:
		t = i * i
		if t == x: return True
		elif t > x: break
		i += 1
	return False

def time_it(func):
	def _(*args, **kwargs):
		import time
		print(str(func).split()[1], 'is running')
		t = time.perf_counter(); r = func(*args, **kwargs)
		t = time.perf_counter() - t
		print(round(t, 3), 'seconds') # print('%.4f' % t)
		return r
	return _

@time_it
def test_1():
	for i in range(1 << 21):
		a = is_square_1(99980001); b = is_square_1(99980000)
	print(a, b)

@time_it
def test_2():
	for i in range(1 << 21):
		a = is_square_2(99980001); b = is_square_2(99980000)
	print(a, b)

@time_it
def test_3():
	for i in range(1 << 21):
		a = is_square_3(99980001); b = is_square_3(99980000)
	print(a, b)

#print(is_square_3(99980001), is_square_3(99980000)); quit()
test_3(); test_2(); test_1()
test_1(); test_2(); test_3()

def fn(n):
	n = 1 << n; mask = n - 1
	s = set()
	for i in range(1<<16): s.add((i*i) & mask)
	l = len(s)
	print('%7d' % l, '%7d' % n, float(l) / n)

for i in range(2, 21): fn(i)
'''
不在表里，一定不是。在表里，不一定是。
      2       4 0.5
      3       8 0.375
      4      16 0.25
      7      32 0.21875
     12      64 0.1875
     23     128 0.1796875
     44     256 0.171875
     87     512 0.169921875
    172    1024 0.16796875
    343    2048 0.16748046875
    684    4096 0.1669921875
   1367    8192 0.1668701171875
   2732   16384 0.166748046875
   5463   32768 0.166717529296875
  10924   65536 0.16668701171875
  21847  131072 0.16667938232421875
  43692  262144 0.1666717529296875
  54615  524288 0.10416984558105469
  60076 1048576 0.057292938232421875
'''

趣闻轶事：1. 2004年，普林斯顿大学的数学家John Conway和Alexander Soifer 向《美国数学月刊》递交了他们认为有史以来字数最少的数学论文。整篇论文的正文就两个字：n² + 2 can. 2. 好像有人连图都不配，就{a number} is not.一句。can确实比is not短啊。