上の方で PARI/GP を推しましたが、今は pythonで何でもできますね。 そういう時代なんですね。
bpython 等のREPL(対話式インターフェース) で、電卓代わりに使えますし、plotも綺麗ですし。
PARI/GPの利点は「ほん少しだけ起動が早い」「追加パッケージのinstall&import が不要」くらいでしょうか。

import time, gmpy2 # gmpy2 はググって見つけた数論パッケージ
t=time.time()
li= list( filter( lambda n:gmpy2.is_prime(n), [10**1000+k for k in range(2000) ]) )
t=time.time()-t
print(f" length: {len(li)}\n time: {t*(10*3):7.2f} msec\n" )
--->
length: 2 {10^1000+0...10^1000+1999 の範囲に素数は2つだけ}
time: 105.01 msec

PARI/GP で同じように 1000桁バージョンやらせたら数分待っても応答無しなので諦めました。
大半の数が直ぐに素数判定できるんですが、判定が苦手な数がポツポツ混じっていました。
isprime(10^1000 + 453) \\ 例えばコレとか... Wolfram先生や python のはアルゴリズムが異なるようですね。