1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677 |
- #!/usr/bin/ruby
- # Smallest strong pseudoprime to base 2 with n prime factors.
- # https://oeis.org/A180065
- # Known terms:
- # 2047, 15841, 800605, 293609485, 10761055201, 5478598723585, 713808066913201, 90614118359482705, 5993318051893040401
- # New terms found (24 September 2022):
- # a(11) = 24325630440506854886701
- # a(12) = 27146803388402594456683201
- # a(13) = 4365221464536367089854499301
- # a(14) = 2162223198751674481689868383601
- # a(15) = 548097717006566233800428685318301
- #`(
- # PARI/GP program (slow):
- strong_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, p, j, k_exp, congr) = my(list=List()); forprime(q=p, sqrtnint(B\m, j), if(base%q != 0, my(tv=valuation(q-1, 2)); if(tv > k_exp && Mod(base, q)^(((q-1)>>tv)<<k_exp) == congr, 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, k_exp, congr)))), break); v *= q; t *= q)))); list); my(res=f(1, 1, 2, k, 0, 1)); for(v=0, logint(B, 2), res=concat(res, f(1, 1, 2, k, v, -1))); vecsort(Vec(res));
- a(n) = if(n < 2, return()); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=strong_fermat_psp(x, y, n, 2)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
- # PARI/GP program (fast):
- strong_check(p, base, e, r) = my(tv=valuation(p-1, 2)); tv > e && Mod(base, p)^((p-1)>>(tv-e)) == r;
- strong_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k, e, r) = 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 && strong_check(p, base, e, r), 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; strong_check(p, base, e, r) || 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, e, r)); q *= p; Mod(base, q)^z == 1 || break; v *= p))); list); my(res=f(1, 1, 2, k, 0, 1)); for(v=0, logint(B, 2), res=concat(res, f(1, 1, 2, k, v, -1))); vecsort(Set(res));
- a(n) = if(n < 2, return()); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=strong_fermat_psp(x, y, n, 2)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
- )
- func a(n) {
- return nil if (n < 2)
- var x = pn_primorial(n+1)/2
- var y = 2*x
- loop {
- #say "Sieving range: #{[x,y]}"
- #var arr = n.strong_fermat_psp(2, x,y)
- var arr = n.squarefree_strong_fermat_psp(2, x,y)
- if (arr.len >= 1) {
- return arr[0]
- }
- x = y+1
- y = 2*x
- }
- }
- for n in (2..20) {
- say "a(#{n}) = #{a(n)}"
- }
- __END__
- ? a(2)
- %3 = 2047
- ? a(3)
- %4 = 15841
- ? a(4)
- %5 = 800605
- ? a(5)
- %6 = 293609485
- ? a(6)
- %7 = 10761055201
- ? a(7)
- %8 = 5478598723585
- ? a(8)
- %9 = 713808066913201
- ? a(9)
- %10 = 90614118359482705
- ? a(10)
- %11 = 5993318051893040401
- ? a(11)
- %12 = 24325630440506854886701
- ?
|