sidef-scripts/Math/squarefree_strong_fermat_pseudoprimes_in_range.sf

#!/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]