sidef-scripts/Math/fermat_overpseudoprimes_in_range.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 04 September 2022
# https://github.com/trizen

# Generate all the k-omega Fermat overpseudoprimes 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

func inverse_znorder_primes(base, lambda) is cached {
    base**lambda - 1 -> prime_divisors
}

func iterate_over_primes(x,y,base,lambda,callback) {

    if ((lambda > 1) && (lambda <= 100)) {
        for q in (inverse_znorder_primes(base, lambda)) {

            next  if (q < x)
            break if (q > y)

            #znorder(base, q) == lambda || next

            callback(q)
        }
        return nil
    }

    if (lambda > 1) {

        var u = lambda*idiv_ceil(x-1, lambda)

        return nil if (u > y)

        for w in (range(u, y, lambda)) {
            with (w+1) { |p|
                callback(p) if p.is_prime
            }
        }
        return nil
    }

    each_prime(x, y, callback)
}

func fermat_overpseudoprimes_in_range(A, B, k, base, callback) {

    A = max(k.pn_primorial, A)

    func (m, lambda, lo, j) {

        var hi = idiv(B,m).iroot(j)

        if (lo > hi) {
            return nil
        }

        iterate_over_primes(lo, hi, base, lambda, {|p|

            next if p.divides(base)

            for (var (q,v) = (p, m*p); v <= B; (q,v) = (q*p, v*p)) {

                var L = znorder(base, q)
                if (lambda > 1) {
                    lambda == L || break
                }
                v.is_coprime(L) || break

                if (j == 1) {
                    v >= A || next
                    k.is_one && v.is_prime && next
                    L `divides` v-1 || next
                    callback(v)
                    next
                }

                __FUNC__(v, L, p+1, j-1)
            }
        })
    }(1, 1, 2, k)

    return callback
}

# Generate all the 3-omega Fermat overpseudoprimes to base 2 in range [1194649, 14325843000]

var k    = 3
var base = 2
var from = 1194649
var upto = 14325843000

say gather { fermat_overpseudoprimes_in_range(from, upto, k, base, { take(_) }) }.sort

__END__
[13421773, 464955857, 536870911, 1220114377, 1541955409, 2454285751, 3435973837, 5256967999, 5726579371, 7030714813, 8493511669, 8538455017, 8788016089, 10545166433, 13893138041]