sidef-scripts/Math/squarefree_fermat_overpseudoprimes_in_range.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 29 August 2022
# https://github.com/trizen

# Generate all the squarefree Fermat overpseudoprimes 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 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 squarefree_fermat_overpseudoprimes_in_range(a, b, k, base, callback) {

    a = max(k.pn_primorial, a)

    func (m, L, lo, k) {

        var hi = idiv(b,m).iroot(k)

        return nil if (lo > hi)

        if (k == 1) {

            lo = max(lo, idiv_ceil(a, m))

            return nil if (lo > hi)

            iterate_over_primes(lo, hi, base, L, {|p|
                if (powmod(base, L, p) == 1) {
                    with (m*p - 1) {|t|
                        if ((L `divides` t) && (L == znorder(base, p))) {
                            callback(t+1)
                        }
                    }
                }
            })

            return nil
        }

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

            p.divides(base) && next

            var z = znorder(base, p)
            if (L > 1) {
                z == L || next
            }
            m.is_coprime(z) || next

            __FUNC__(m*p, z, p.next_prime, k-1)
        })
    }(1, 1, 2, k)

    return callback
}

# High-level implementation (unoptimized)
func squarefree_fermat_overpseudoprimes_in_range_unoptimized(A, B, k, base, callback) {

    func F(m, lambda, p, k) {
        if (k == 1) {
            each_prime(max(p, ceil(A/m)), floor(B/m), {|q|
                if ((powmod(base, lambda, q) == 1) && (znorder(base, q) == lambda) && (lambda `divides` (m*q - 1))) {
                    callback(m*q)
                }
            })
        }
        elsif (k >= 2) {
            for q in (primes(p, (B/m)**(1/k))) {
                if (!q.divides(base)) {
                    var L = znorder(base, q)
                    if ((lambda==L || lambda==1) && (gcd(L, m) == 1)) {
                        F(m*q, L, q.next_prime, k-1)
                    }
                }
            }
        }
    }

    F(1, 1, 2, k)
}

# Generate all the squarefree Fermat overpseudoprimes to base 2 with 3 prime factors in the range [1194649, 14325843000]

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

say gather { squarefree_fermat_overpseudoprimes_in_range(from, upto, k, base, { take(_) }) }.sort
say gather { squarefree_fermat_overpseudoprimes_in_range_unoptimized(from, upto, k, base, { take(_) }) }.sort if false

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