sidef-scripts/Math/squarefree_lucas_U_pseudoprimes_in_range.sf

#!/usr/bin/ruby

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

# Generate all the squarefree Lucas pseudoprimes to the U_n(P,Q) sequence with n prime factors in a given range [a,b]. (not in sorted order)

# See also:
#   https://en.wikipedia.org/wiki/Almost_prime
#   https://en.wikipedia.org/wiki/Lucas_sequence
#   https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html

func lucas_znorder(P,Q,D,n) {
    n - kronecker(D, n) -> divisors.first_by {|d| lucasUmod(P, Q, d, n) == 0 }
}

func squarefree_lucas_U_pseudoprimes_in_range(a, b, k, P, Q, callback) {

    a = max(k.pn_primorial, a)

    var D = (P*P - 4*Q)

    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))
            lo > hi && return nil

            for j in (1, -1) {
                var t = mulmod(m.invmod(L), j, L)

                t > hi && next
                t += L*idiv_ceil(lo - t, L) if (t < lo)
                t > hi && next

                for p in (range(t, hi, L)) {
                    p.is_prime || next
                    with (m*p) {|n|
                        with (n - kronecker(D, n)) {|w|
                            if ((L `divides` w) && (lucas_znorder(P, Q, D, p) `divides` w)) {
                                callback(n)
                            }
                        }
                    }
                }
            }

            return nil
        }

        each_prime(lo, hi, {|p|

            p.divides(D) && next

            var z = lucas_znorder(P, Q, D, p)
            m.is_coprime(z) || next

            __FUNC__(m*p, lcm(L, z), p+1, k-1)
        })
    }(1, 1, 2, k)

    return callback
}

# Generate all the squarefree Fibonacci pseudoprimes in the range [1, 15251]

var from = 1
var upto = 15251
var P    = 1
var Q    = -1
var arr  = []

for k in (2..100) {
    break if k.pn_primorial>upto
    squarefree_lucas_U_pseudoprimes_in_range(from, upto, k, P, Q, { arr << _ })
}

say arr.sort

__END__
[323, 377, 1891, 3827, 4181, 5777, 6601, 6721, 8149, 10877, 11663, 13201, 13981, 15251]