Pythonで素数A＋素数B＝1000000となる素数の組をすべて探してみた

はじめに

表題にあるように素数A＋素数B＝1000000となる素数の組をすべて探します。例えば、素数Aに17、素数Bに999987を取ると1000000となります。簡単に書いたので、読みにくいと思います。特に変数名...。ご了承ください。

実装

まずは、1000000までの素数をすべて探してリストに保存していきます。ある数ｂをこれまでに見つかった素数で割っていき、√b以下の素数すべてで割り切れないならば、素数リストに追加していくというアルゴリズムです。ついでに素数の個数も数えました。

# !/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

list_prime_N=[2,3]
a=1
b=3	
k=1000000       

while b<=k:
		c=list_prime_N[a]
		extra=b%c
		if extra==0:
			b+=2
			a=0
		else:
			a+=1
			d=np.sqrt(b)
			if list_prime_N[a]>d:
				list_prime_N=list_prime_N+[b]
				b+=2
				a=0

print(list_prime_N)
print(len(list_prime_N))

次に足して1000000になる数を探していきます。list_prime_Nには、1000000までの素数が昇順に入っているので、まずは１番大きい素数と小さい素数を足してみます。1000000より大きければ、１番大きな素数を2番めに大きい素数に変え、1000000より小さければ、1番目に小さな素数を2番めに小さな素数に変えます。これをくりかえしていきます。つまり1000000を超えれば、大きい方から小さい方へ、1000000を下回れば小さい方から大きい方へ動かしていきます。これは、コードを見たほうがわかりやすいと思います。昇順に動いてきた方と降順に動いてきたほうが入れ替えるまで繰り返します。ついでに組みの数も数えました。

A=0
B=-1
cnt=0
P=list_prime_N[A]
Q=list_prime_N[B]
pp=[[]]
while P<Q:
	P=list_prime_N[A]
	Q=list_prime_N[B]
	if P+Q<k:
		A+=1
	elif P+Q==k:
		pp+=[[P,Q]]
		cnt+=1
		A+=1
		B-=1
	elif P+Q>k:
		B=B-1
		
print(pp)
print(cnt)

---

このプログラムによると1000000までの素数は78498個、組の総数は5402組のようです。このプログラムのk（1000000）を変えていくとゴールドバッハの予想の検証になります。また、このプログラムはシングルコアで計算していくので（しかもpython）速度はあまり早くありません。
間違えがあったら教えて下さい。