1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 07 November 2022
- # https://github.com/trizen
- # Generate all the squarefree strong Fermat pseudoprimes to a given base with n prime factors in a given range [a,b]. (not in sorted order)
- # See also:
- # https://en.wikipedia.org/wiki/Almost_prime
- # https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
- func squarefree_strong_fermat_pseudoprimes_in_range(A, B, k, base, callback) {
- A = max(k.pn_primorial, A)
- var generator = func (m, L, lo, k, k_exp, congr) {
- var hi = idiv(B,m).iroot(k)
- return nil if (lo > hi)
- if (k == 1) {
- lo = max(lo, idiv_ceil(A, m))
- lo > hi && 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 || next
- with (m*p) {|n|
- if (znorder(base, p) `divides` n.dec) {
- var v = p.dec.valuation(2)
- if (v > k_exp && powmod(base, p.dec>>(v-k_exp), p).is_congruent(congr, p)) {
- callback(n)
- }
- }
- }
- }
- return nil
- }
- each_prime(lo, hi, {|p|
- p.divides(base) && next
- var val2 = p.dec.valuation(2)
- val2 > k_exp || next
- powmod(base, p.dec>>(val2-k_exp), p).is_congruent(congr, p) || next
- var z = znorder(base, p)
- m.is_coprime(z) || next
- __FUNC__(m*p, lcm(L, z), p+1, k-1, k_exp, congr)
- })
- }
- # Case where 2^d == 1 (mod p), where d is the odd part of p-1.
- generator(1, 1, 2, k, 0, 1)
- # Cases where 2^(d * 2^v) == -1 (mod p), for some v >= 0.
- for v in (0..B.ilog2) {
- generator(1, 1, 2, k, v, -1)
- }
- return callback
- }
- # Generate all the squarefree strong Fermat pseudoprimes to base 5 with 4 prime factors in the range [1, 10^8]
- var k = 4
- var base = 5
- var from = 1
- var upto = 1e8
- say gather { squarefree_strong_fermat_pseudoprimes_in_range(from, upto, k, base, { take(_) }) }.sort
- __END__
- [4382191, 7267051, 9694351, 11921001, 18464761, 28351081, 33333091, 33567451, 38574199, 42151351, 48321001, 71107681, 80717131, 83420101]
|