无心人 发表于 2020-12-1 19:00
我感觉模乘不需要除法
最后取模你准备怎么搞？
我已经把取模的数量精简到一个了
mulmod=|mut x,mut y|{
let =;
x=lx*hy+ly*hx;
y=lx*ly+((x&9223372036854775807)<<63);x>>=63;//y=low_mid, x=high_mid<<126
//lx*ly+y+x<<126+hx*hy<<126
x+=hx*hy;
//lx*ly+y+x<<126
hx=x&2147483647;
x>>=31;
//(y+hx*m126+x*m157)%module
y+=hx*m126+x*m157;
//y-(((y>>64)as u64)/m64) as u128*module
y%module
//lx.max(ly) < 1<<63 ==> lx*ly < 1<<126
// y < 9223372036854775807<<63 < 1<<126
//hx <2147483647*module< 1<<31+94 = 1<<126
// x <(module*module>>(63*2+31))*module < 1<<(94*2-63*2-31+94)=1<<125
//Thus, lx*ly+y+hx*m126+x*m157 < 1<<128, safe to add them together and executing % ops.
};
刚查了gmp的源代码
mpz_millerrabin(mpz_t n, int d);
d是测试次数，而不是测试基
所以不合适，需要自己写
说下只用加减乘的模乘
u=u1u0
v=v1v0
w=w1w0
r0=r3r2r1r0=u*v
r1=r0-(r3r2)*w<r0
这样设计是不是可以
#include<gmp.h>
#include<stdio.h>
#include<stdlib.h>
#include<stdbool.h>
#define NULL_TERMINATED0
bool fermatTest( const mpz_t n, unsigned long b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
int result = mpz_cmp_ui(r, 1);
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return (result == 0);
}
bool millerrabinTest( const mpz_t n, unsigned long b )
{
mpz_t ns1, t, x, b1, r;
unsigned int s;
bool result;
mpz_inits(ns1, t, x, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);//b1=b
mpz_sub_ui(ns1, n, 1);//ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);//ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (int i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
mpz_clears(ns1, t, x, b1, r, NULL_TERMINATED);
return result;
}
void testFunc( void )
{
mpz_t n;
int d;
mpz_init(n);
mpz_set_ui(n, 341); d = 2;
gmp_printf("fermat test(%Zd, %d) = %d\n", n, d, fermatTest(n, d));
gmp_printf("miller rabin test(%Zd, %d) = %d\n", n, d, millerrabinTest(n, d));
mpz_set_ui(n, 8321); d = 2;
gmp_printf("fermat test(%Zd, %d) = %d\n", n, d, fermatTest(n, d));
gmp_printf("miller rabin test(%Zd, %d) = %d\n", n, d, millerrabinTest(n, d));
mpz_set_ui(n, 65541); d = 2;
gmp_printf("fermat test(%Zd, %d) = %d\n", n, d, fermatTest(n, d));
gmp_printf("miller rabin test(%Zd, %d) = %d\n", n, d, millerrabinTest(n, d));
mpz_set_ui(n, 65539); d = 3;
gmp_printf("fermat test(%Zd, %d) = %d\n", n, d, fermatTest(n, d));
gmp_printf("miller rabin test(%Zd, %d) = %d\n", n, d, millerrabinTest(n, d));
mpz_set_ui(n, 65539); d = 3;
mpz_pow_ui(n, n, 8);
mpz_add_ui(n, n, 78);
gmp_printf("fermat test(%Zd, %d) = %d\n", n, d, fermatTest(n, d));
gmp_printf("miller rabin test(%Zd, %d) = %d\n", n, d, millerrabinTest(n, d));
mpz_clear(n);
}
int main( void )
{
testFunc();
return 0;
}
fermat test(341, 2) = 1
miller rabin test(341, 2) = 0
fermat test(8321, 2) = 1
miller rabin test(8321, 2) = 1
fermat test(65541, 2) = 0
miller rabin test(65541, 2) = 0
fermat test(65539, 3) = 1
miller rabin test(65539, 3) = 1
fermat test(340407002012868253357199029315156056559, 3) = 1
miller rabin test(340407002012868253357199029315156056559, 3) = 1
无心人 发表于 2020-12-2 10:36
341 2 1
65541 2 0
65539 3 1
你可能需要这个东西：
gmp-6.2.1/mini-gmp/mini-gmp.c
static int
gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
const mpz_t q, mp_bitcnt_t k)
{
assert (k > 0);
/* Caller must initialize y to the base. */
mpz_powm (y, y, q, n);
if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
return 1;
while (--k > 0)
{
mpz_powm_ui (y, y, 2, n);
if (mpz_cmp (y, nm1) == 0)
return 1;
/* y == 1 means that the previous y was a non-trivial square root
of 1 (mod n). y == 0 means that n is a power of the base.
In either case, n is not prime. */
if (mpz_cmp_ui (y, 1) <= 0)
return 0;
}
return 0;
}
2.5G主频的处理器，gmp测试10^8的时间是142.8秒
3.2G主频的处理器，同样逻辑的Julia的时间是236.5秒
另外发现，gmp的gmp_printf在windows下会崩，找不到崩的原因，所以只好少用~
#include<gmp.h>
#include<stdbool.h>
#include<stdio.h>
#include<stdlib.h>
#define NULL_TERMINATED 0
bool fermatTest( const mpz_t n, unsigned long b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
bool result = mpz_cmp_ui(r, 1) == 0;
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return result;
}
bool millerrabinTest( const mpz_t n, unsigned long b )
{
mpz_t ns1, t, b1, r;
unsigned long s;
bool result;
mpz_inits(ns1, t, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);// b1=b
mpz_sub_ui(ns1, n, 1);// ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);// ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (unsigned long i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
mpz_clears(ns1, t, b1, r, NULL_TERMINATED);
return result;
}
const unsigned long step = 100000000;
// 实际测试值=base + index*step +
void init( mpz_t base, unsigned long index )
{
mpz_t tmp;
mpz_set_ui(base, 10);
mpz_pow_ui(base, base, 28);
mpz_init_set_ui(tmp, step);
mpz_mul_ui(tmp, tmp, index);
mpz_add(base, base, tmp);
mpz_clear(tmp);
}
bool quickTest( mpz_t n )
{
unsigned long b[] = {3, 5, 7, 11, 13};
for (int i = 0; i < sizeof(b)/sizeof(unsigned long); i++)
if (!millerrabinTest(n, b))
return false;
return true;
}
void testFunc( mpz_t base )
{
mpz_t n;
for (unsigned long i = 1; i < step; i += 2)
{
mpz_init(n);
mpz_add_ui(n, base, i);
if (fermatTest(n, 2))
if (millerrabinTest(n, 2))
{
if (!quickTest(n))
gmp_printf("%Zd : SPSP(2)\n", n);
}
else
gmp_printf("%Zd : PSP(2)\n", n);
mpz_clear(n);
}
}
int main( void )
{
mpz_t base;
mpz_init(base);
init( base, 0 );
testFunc( base );
mpz_clear(base);
return 0;
}
using Primes
global function fermatProbablePrimeTest(n, b)
a = powermod(b, n-1, n)
return (a == 1)
end
global function millerRabinProbablePrimeTest( n, b )
k = trailing_zeros(n-1)
m = (n-1) >> k
a = powermod(b, m, n)
if (a == 1) || (a == (n-1))
return true
else
for i in 1:k-1
a = powermod(a, 2, n)
if a == n - 1
return true
end
end
end
return false
end
global function quickPrimeTest(n)
for b in
if !millerRabinProbablePrimeTest(n, b)
return false
end
end
return true
end
global b28 = big(10)^28
println("==========")
@time begin
for i in range(1, step = 2, stop = 99999999)
test = b28 + i
if fermatProbablePrimeTest(test, 2)
if millerRabinProbablePrimeTest(test, 2)
if !quickPrimeTest(test)
println("$test is SPSP(2)")
end
else
println("$test is PSP(2)")
end
end
end
end
println("==========")
无心人 发表于 2020-12-4 09:38
2.5G主频的处理器，gmp测试10^8的时间是142.8秒
3.2G主频的处理器，同样逻辑的Julia的时间是236.5秒
./test
free(): double free detected in tcache 2
已放弃 (核心已转储)
看上去像是93和95行把同一个n clear了两次的结果
gcc test.c -march=native -O3 -pipe -flto -lgmp -o test && time ./test
real 1m20.758s
user 1m20.706s
sys 0m0.000s
在我这里测试的结果只是比Rust略快，而快的原因可以解释成，这里用的是quickTest而非ispseudoprime
感觉你大概重复造了跟轮子……
装一个Rust不好吗：）
做了筛法的一个测试，1亿的数字需要测试fermat（2）的数字从5000万降低到大概1200万~
还没确认程序正确，代码可能会修改
粗略测试，20核并发，1亿数字6秒左右
#include<gmp.h>
#include<omp.h>
#include<stdbool.h>
#include<stdio.h>
#include<stdlib.h>
#define NULL_TERMINATED 0
bool fermatTest( const mpz_t n, unsigned long b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
bool result = mpz_cmp_ui(r, 1) == 0;
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return result;
}
bool millerrabinTest( const mpz_t n, unsigned long b )
{
mpz_t ns1, t, b1, r;
unsigned long s;
bool result;
//gmp_printf("MR(%Zd, %d)...\n", n, b);
mpz_inits(ns1, t, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);// b1=b
mpz_sub_ui(ns1, n, 1);// ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);// ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (unsigned long i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
//gmp_printf("MR(%Zd, %d) = %d\n", n, b, result);
mpz_clears(ns1, t, b1, r, NULL_TERMINATED);
return result;
}
const unsigned long step = 100000000;
char*buff;
unsigned long cnt = 0;
unsigned long small;
unsigned long p100[] = { 2,3,5,7, 11,
13, 17, 19, 23, 29,
31, 37, 41, 43, 47,
53, 59, 61, 67, 71,
73, 79, 83, 89, 97};
//p < 10000
unsigned long ord(unsigned long p)
{
unsigned long factor;
unsigned long s = 0;
unsigned long n = p - 1;
mpz_t b, tmp, mp;
mpz_init_set_ui(b, 2);
mpz_init_set_ui(mp, p);
mpz_init(tmp);
for (int i = 0; i < 25; i ++)
{
if (n%p100 == 0)
do
{
factor = p100;
n /= p100;
} while (n%p100 == 0);
}
n = p - 1;
for (int i = 0; i < s; i++)
{
mpz_powm_ui(tmp, b, n / factor, mp);
if (mpz_cmp_ui(tmp, 1) == 0)
n /= factor;
}
if (n%2==1)
n*=2;
mpz_clears(b, tmp, mp, NULL_TERMINATED);
return n;
}
// 实际测试值=base + index*step +
void init( mpz_t base, unsigned long index )
{
mpz_t tmp, tmp1;
mpz_set_ui(base, 10);
mpz_pow_ui(base, base, 28);
mpz_init_set_ui(tmp, step);
mpz_mul_ui(tmp, tmp, index);
mpz_add(base, base, tmp);
buff = (char *)malloc(step);
for (long i = 0; i < step; i += 2)
buff = 0;
for (long i = 1; i < step; i += 2)
buff = 1;
unsigned long i = 0;
mpz_init_set_ui(tmp1, 3);
while (i < 1227)
{
mpz_nextprime(tmp, tmp1);
mpz_set(tmp1, tmp);
small = mpz_get_ui(tmp);
small = ord(small);
i++;
}
mpz_clears(tmp, tmp1, NULL_TERMINATED);
}
bool quickTest( mpz_t n )
{
//gmp_printf("%Zd quick test...\n", n); fflush(stdout);
unsigned long b[] = {3, 5, 7, 11, 13};
for (int i = 0; i < sizeof(b)/sizeof(unsigned long); i++)
if (!millerrabinTest(n, b))
return false;
return true;
}
void check(mpz_t base, long i)
{
mpz_t n;
mpz_init(n);
mpz_add_ui(n, base, i);
if (fermatTest(n, 2))
if (millerrabinTest(n, 2))
{
//gmp_printf("quick test.....\n");
if (!quickTest(n))
#pragma omp critical
{
gmp_printf("%Zd : SPSP(2)\n", n);
}
}
else
#pragma omp critical
{
gmp_printf("%Zd : PSP(2)\n", n);
}
mpz_clear(n);
#pragma omp atomic
cnt++;
}
void testFunc( mpz_t base )
{
mpz_t q, r, q1;
mpz_inits(q, r, q1, NULL_TERMINATED);
gmp_printf("start: %Zd\n", base);
for (long i = 0; i < 1227; i ++)
{
mpz_fdiv_qr_ui(q, r, base, small);
mpz_mod_ui(q1, q, small);
unsigned long q2 = mpz_get_ui(q1);
unsigned long r1 = mpz_get_ui(r);
r1 = small - r1;
q2 = (q2+1)%small;
if (r1%2==0)
{
r1 += small;
q2 = (q2+1)%small;
}
while (r1 < step)
{
if (q2 != 1)
buff = 0;
r1 += 2*small;
q2 = (q2+2)%small;
}
}
#pragma omp parallel for schedule(static,20)
for (long i = 1; i < step; i += 2)
if (buff == 1)
check(base, i);
mpz_clears(q, r, q1, NULL_TERMINATED);
}
void clear( void )
{
free(buff);
}
int main( void )
{
mpz_t base;
mpz_init(base);
init( base, 0 );
testFunc( base );
printf("%d fermat test\n", cnt);
mpz_clear(base);
clear();
return 0;
}
100万范围的测试，50万个奇数
2.5G主频处理器，单线程
一次Fermat测试，6次Miller Robin测试
C 1.18秒， Julia 2.65秒， go 6.46秒
C
#include<gmp.h>
#include<stdbool.h>
#include<stdio.h>
#include<stdlib.h>
#define NULL_TERMINATED 0
bool fermatTest( const mpz_t n, unsigned long b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
bool result = mpz_cmp_ui(r, 1) == 0;
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return result;
}
bool millerrabinTest( const mpz_t n, unsigned long b )
{
mpz_t ns1, t, b1, r;
unsigned long s;
bool result;
mpz_inits(ns1, t, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);// b1=b
mpz_sub_ui(ns1, n, 1);// ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);// ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (unsigned long i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
mpz_clears(ns1, t, b1, r, NULL_TERMINATED);
return result;
}
const long step = 1000000;
void init( mpz_t base, int idx)
{
mpz_t tmp;
mpz_set_ui(base, 10);
mpz_pow_ui(base, base, 28);
mpz_init_set_ui(tmp, step);
mpz_mul_ui(tmp, tmp, idx);
mpz_add(base, base, tmp);
mpz_clears(tmp, NULL_TERMINATED);
}
bool quickTest( mpz_t n )
{
//gmp_printf("%Zd quick test...\n", n); fflush(stdout);
unsigned long b[] = {3, 5, 7, 11, 13};
for (int i = 0; i < sizeof(b)/sizeof(unsigned long); i++)
if (!millerrabinTest(n, b))
return false;
return true;
}
void check(mpz_t base, long i)
{
mpz_t n;
mpz_init(n);
mpz_add_ui(n, base, i);
if (fermatTest(n, 2))
if (millerrabinTest(n, 2))
{
if (!quickTest(n))
{
gmp_printf("%Zd : SPSP(2)\n", n);
}
}
else
{
gmp_printf("%Zd : PSP(2)\n", n);
}
mpz_clear(n);
}
void testFunc( mpz_t base )
{
mpz_t q, r, q1;
mpz_inits(q, r, q1, NULL_TERMINATED);
for (long i = 1; i < step; i += 2)
check(base, i);
mpz_clears(q, r, q1, NULL_TERMINATED);
}
int main( int argc, char * argv[] )
{
mpz_t base;
mpz_init(base);
init( base, 0);
gmp_printf("%Zd start...\n", base);
testFunc( base );
mpz_clear(base);
return 0;
}
Julia
using Primes
global function fermatTest(n, b)
a = powermod(b, n-1, n)
return (a == 1)
end
global function millerRabinTest( n, b )
k = trailing_zeros(n-1)
m = (n-1) >> k
a = powermod(b, m, n)
if (a == 1) || (a == (n-1))
return true
else
for i in 1:k-1
a = powermod(a, 2, n)
if a == n - 1
return true
end
end
end
return false
end
global function quickTest(n)
for b in
if !millerRabinTest(n, b)
return false
end
end
return true
end
global b28 = big(10)^28
@time begin
for i in range(1, step = 2, stop = 999999)
test = b28 + i
if fermatTest(test, 2)
if millerRabinTest(test, 2)
if !quickTest(test)
println("$test is SPSP(2)")
end
else
println("$test is PSP(2)")
end
end
end
end
Go
package main
import (
"fmt"
"math/big"
"time"
)
func fermatTest(n *big.Int, b int64) bool {
ns1 := big.NewInt(0)
ns1.Set(n)
one := big.NewInt(1)
b1 := big.NewInt(b)
ns1.Sub(ns1, one)
ns1.Exp(b1, ns1, n)
return ns1.Cmp(one) == 0
}
func millerRabinTest(n *big.Int, b *big.Int) bool {
s := int(0)
ns1 := big.NewInt(0)
ns1.Set(n)
one := big.NewInt(1)
two := big.NewInt(2)
ns1.Sub(ns1, one)
t := big.NewInt(0)
t.Set(ns1)
for t.Bit(0) == 0 {
s++
t.Rsh(t, 1)
}
t.Exp(b, t, n)
if t.Cmp(one) == 0 || t.Cmp(ns1) == 0 {
return true
}
for i := 1; i < s; i++ {
t.Exp(t, two, n)
if t.Cmp(ns1) == 0 {
return true
}
}
return false
}
func quickTest(n *big.Int) bool {
for _, t := range []int64{3, 5, 7, 11, 13} {
b := big.NewInt(t)
if !millerRabinTest(n, b) {
return false
}
}
return true
}
func main() {
t0 := time.Now()
b := big.NewInt(10)
n28 := big.NewInt(28)
two := big.NewInt(2)
b.Exp(b, n28, nil)
for i := 1; i < 1000000; i += 2 {
t := big.NewInt(int64(i))
t.Add(b, t)
if fermatTest(t, 2) {
if millerRabinTest(t, two) {
if !quickTest(t) {
fmt.Println(t.String(), " is SPSP(2)")
}
} else {
fmt.Println(t.String(), " is PSP(2)")
}
}
}
fmt.Println(time.Since(t0).String())
}
本帖最后由 .·.·. 于 2020-12-8 21:40 编辑
rust 优化版from=10000000000000000000000000001,len=250
10000000000000000000100000001.1552676,0,34.23762839s搜索前一亿个数字，不开并行，34.23762839s（得到1552676个素数，与之前的程序一致）使用nightly版本的rust，代码如下：#!
use rug::{Assign, Integer};
const SIEVE_SIZE:usize=200000;
const PRIME_SIZE:usize=include!("final.txt").len();
const PRIME_MAP:[;PRIME_SIZE]=include!("final.txt");
fn main() {
const CONSTARR:=;
let now=std::time::Instant::now();
let arguments: Vec<String> = std::env::args().collect();
let begin=if let Some(x)=arguments.get(2){x.parse::<u128>().unwrap_or(0)}else{0};
let len=if let Some(x)=arguments.get(1){x.parse::<u128>().unwrap_or(250)}else{250};
let start=10000000000000000000000000000u128+begin*len;
let mut am1=Integer::from(start);
let mut a=Integer::from(start+1);
println!("from={},len={}",a,len);
let mut b=Integer::from(2);
let mut am1d=Integer::from(0);
let mut prime=0u64;
let mut psp=0u64;
let mut sieve=;
let mut prime_curr=;
let mut prime_order=;
for i in 0..PRIME_SIZE{
let t=(start%PRIME_MAP as u128) as i32;
let next= if t==0{
start+PRIME_MAP as u128
}else{
let mut tmp=start-(t as u128)+PRIME_MAP as u128;
if tmp%2==0{tmp+=PRIME_MAP as u128}
tmp
};
let count=next/PRIME_MAP as u128;
prime_curr=((next-start)>>1) as i32;
prime_order=PRIME_MAP-((count>>1)%PRIME_MAP as u128) as i32
}
for _i in 0..len{
sieve.fill(true);
for (i,&) in PRIME_MAP.iter().enumerate(){
let mut cur=prime_curr;
let mut ord=prime_order;
while cur<SIEVE_SIZE as i32{
if ord==0{
ord=order
}else{
sieve=false
}
cur+=prime;
ord-=1
}
prime_curr=cur-SIEVE_SIZE as i32;
prime_order=ord
}
for test in &sieve{
if *test{
let aa=am1.find_one(0).unwrap().min(6);
am1d.assign(&am1>>aa);
if {b.assign(CONSTARR);let _=b.pow_mod_mut(&am1d,&a);b==1u8}{//Fermat
psp+=1;
if {b.assign(2u8);//Miller-Rabin
let aa=a.find_one(0).unwrap();
am1d.assign(&am1>>aa);
let _=b.pow_mod_mut(&am1d,&a);
(b==1 || b==am1 )||{
let mut f=false;
for _ in 0..aa{
let _=b.pow_mod_mut(&Integer::from(2),&a);
if b==am1{f=true;break}
}
f
}
}{
if a.is_probably_prime(25)==rug::integer::IsPrime::No{
println!("{} is SPSP(2)",a)
}else{prime+=1}
}else{println!("{} is PSP(2)",a)}
}
}
a+=2;
am1+=2;
}
}
println!("{}.{},{},{:?}",a,prime,prime-psp,now.elapsed());
}其中final.txt是大于3小于100000的全部9590个素数，以及其order或其order的一半（如果order能被2整除）
可以用gp+python生成这个order列表（gp生成python排版，不换行的话对眼睛和文本编辑器都不怎么友好）
生成用代码如下'''
#generate by GP:
echo "print(vector(primepi(100000)-2,i,my(a=prime(i+2));b=znorder(Mod(2,a));))" | gp -fq > primev2.txt
#processing by python:
with open('primesv2.txt') as f:g=eval(f.read())
#finally print result:
with open('final.txt','w') as f:
_=f.write('[')
t=0
for i in g:
_=f.write(str(i)+',')
t+=1
if t==10:
t=0
_=f.write('\n')
_=f.write(']')
'''成品：[,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
,,,,,,,,,,
]update：如果增加平方项的筛选（删除全部平方数），可以稍快一点
只要不是1093跟3511这两个数字，其他（适合初筛的）p计算mod p^2的阶的时候，阶一定能被p整除，也就是，如果$2^{p^2q-1}=1(mod p^2)$，一定有$p|p^2q-1$也就是$p|-1$，这是不可能的。
所以凡是平方数都可以直接删掉
顺便改了一个（似乎没影响的）BUG
我写的sieve_size=200000是sieve的长度，等效的size其实是400000，现在将sieve_size_div_2设置为真正的sieve长度，而sieve_size弃用。
这里的版本看上去很难继续优化了。
就这样吧