123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 30 August 2022
- # https://github.com/trizen
- # Generate all the k-omega Fermat pseudoprimes in range [a,b]. (not in sorted order)
- # Definition:
- # k-omega primes are numbers n such that omega(n) = k.
- # See also:
- # https://en.wikipedia.org/wiki/Almost_prime
- # https://en.wikipedia.org/wiki/Prime_omega_function
- # https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
- # 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)));
- # 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)));
- func fermat_pseudoprimes_in_range(A, B, k, base, callback) {
- A = max(k.pn_primorial, A)
- var seen = Hash()
- func (m, L, lo, j) {
- var hi = idiv(B,m).iroot(j)
- lo > hi && return nil
- if (j == 1) {
- if (L == 1) { # optimization
- each_prime(lo, hi, {|p|
- p.divides(base) && next
- for (var v = (m == 1 ? p*p : m*p); v <= B; v *= p) {
- v >= A || next
- v.is_psp(base) || break
- callback(v) if !(seen{v} := 0 ++)
- }
- })
- return nil
- }
- var t = m.invmod(L)
- t > hi && return nil
- t += L*idiv_ceil(lo - t, L) if (t < lo)
- t > hi && return nil
- for p in (range(t, hi, L)) {
- p.is_prime_power || next
- p.is_coprime(m) || next
- p.is_coprime(base) || next
- var v = m*p
- v >= A || next
- if (znorder(base, p) `divides` v.dec) {
- callback(v) if !(seen{v} := 0 ++)
- }
- }
- return nil
- }
- each_prime(lo, hi, {|p|
- p.divides(base) && next
- var z = znorder(base, p)
- m.is_coprime(z) || next
- var v = m*p
- var q = p
- while (v <= B) {
- __FUNC__(v, lcm(L, z), p+1, j-1)
- q *= p
- powmod(base, z, q) == 1 || break
- v *= p
- }
- })
- }(1, 1, 2, k)
- return callback
- }
- # Generate all the Fermat pseudoprimes to base 3 in range [1, 10^4]
- var from = 1
- var upto = 1e4
- var base = 3
- var arr = gather {
- for k in (1..100) {
- break if (k.pn_primorial > upto)
- fermat_pseudoprimes_in_range(from, upto, k, base, { take(_) })
- }
- }
- say arr.sort
- # Run some tests
- if (false) { # true to run some tests
- for k in (1..5) {
- say "Testing: k = #{k}"
- var lo = k.pn_primorial
- var hi = lo*1000
- for base in (2..100) {
- var this = k.fermat_psp(base, lo, hi)
- #var this = k.omega_primes(lo,hi).grep{.is_psp(base) && .is_composite }
- var that = gather {
- fermat_pseudoprimes_in_range(lo, hi, k, base, func (n) { take(n) })
- }.sort
- this == that ||
- die "Error for k = #{k} and base = #{base} with hi = #{hi}\n#{this} != #{that}";
- }
- }
- }
- __END__
- [91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601, 7381, 8401, 8911]
|