-
Notifications
You must be signed in to change notification settings - Fork 0
/
prog.sf
53 lines (39 loc) · 2.53 KB
/
prog.sf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#!/usr/bin/ruby
# a(n) = smallest pseudoprime to base 2 with n prime factors.
# https://oeis.org/A007011
# Known terms:
# 341, 561, 11305, 825265, 45593065, 370851481, 38504389105, 7550611589521, 277960972890601, 32918038719446881, 1730865304568301265, 606395069520916762801, 59989606772480422038001, 6149883077429715389052001, 540513705778955131306570201, 35237869211718889547310642241
# New terms found:
# a(18) = 13259400431578770557375638157521
# a(19) = 580827911915963785382253469211401
# a(20) = 292408776547176479576081475390697801
# a(21) = 39498823114155235923831808284152901601
# a(22) = 3284710806953570725820888298298841565601
# a(23) = 327373784481535488655521620744179013043601
# a(24) = 221404014114397213806721960178887462402717201
# a(25) = 43691666165877828056799483424028795272585383601
# a(26) = 13213974925373194568934435211730355813060799098001
# a(27) = 1952204134080476076724242017017925744953021675628161
#`(
# PARI/GP program (slow):
fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, p, j) = my(list=List()); forprime(q=p, sqrtnint(B\m, j), if(base%q != 0, my(v=m*q, t=q, r=nextprime(q+1)); while(v <= B, my(L=lcm(l, znorder(Mod(base, t)))); if(gcd(L, v) == 1, if(j==1, if(v>=A && if(k==1, !isprime(v), 1) && (v-1)%L == 0, listput(list, v)), if(v*r <= B, list=concat(list, f(v, L, r, j-1)))), break); v *= q; t *= q))); list); vecsort(Vec(f(1, 1, 2, k)));
a(n) = if(n < 2, return()); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=fermat_psp(x, y, n, 2)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
# PARI/GP program (fast):
fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, forstep(p=lift(1/Mod(m, l)), hi, l, if(isprimepower(p) && gcd(m*base, p) == 1, my(n=m*p); if(n >= A && (n-1) % znorder(Mod(base, p)) == 0, listput(list, n)))), forprime(p=lo, hi, base%p == 0 && next; my(z=znorder(Mod(base, p))); gcd(m,z) == 1 || next; my(q=p, v=m*p); while(v <= B, list=concat(list, f(v, lcm(l, z), p+1, k-1)); q *= p; Mod(base, q)^z == 1 || break; v *= p))); list); vecsort(Set(f(1, 1, 2, k)));
a(n) = if(n < 2, return()); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=fermat_psp(x, y, n, 2)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
)
func T(k, base=2) {
var x = pn_primorial(k)
var y = 3*x
loop {
var arr = k.fermat_psp(base, x, y)
if (arr) {
return arr[0]
}
x = y+1
y = 3*x
}
}
for n in (2..15) {
say "a(#{n}) = #{T(n)}"
}