sidef-scripts/Math/carmichael_strong_fermat_pseudoprimes_in_range.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 07 November 2022
# https://github.com/trizen

# Generate all the Carmichael numbers with n prime factors in a given range [A,B] that are also strong Fermat pseudoprimes to a given base. (not in sorted order)

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

# PARI/GP program:
#   carmichael_strong_psp(A, B, k, base) = A=max(A, vecprod(primes(k+1))\2); (f(m, l, p, k, k_exp, congr, u=0, v=0) = my(list=List()); if(k==1, forprime(q=u, v, my(t=m*q); if((t-1)%l == 0 && (t-1)%(q-1) == 0, my(tv=valuation(q-1, 2)); if(tv > k_exp && Mod(base, q)^(((q-1)>>tv)<<k_exp) == congr, listput(list, t)))), forprime(q = p, sqrtnint(B\m, k), if(base%q != 0, my(tv=valuation(q-1, 2)); if(tv > k_exp && Mod(base, q)^(((q-1)>>tv)<<k_exp) == congr, my(L=lcm(l, q-1)); if(gcd(L, m) == 1, my(t = m*q, u=ceil(A/t), v=B\t); if(u <= v, my(r=nextprime(q+1)); if(k==2 && r>u, u=r); list=concat(list, f(t, L, r, k-1, k_exp, congr, u, v)))))))); list); my(res=f(1, 1, 3, k, 0, 1)); for(v=0, logint(B, 2), res=concat(res, f(1, 1, 3, k, v, -1))); vecsort(Vec(res));

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

    A = max(k.pn_primorial, A)

    # Largest possisble prime factor for Carmichael numbers <= B
    var max_p = ((1 + isqrt(8*B + 1))>>2)

    var generator = func (m, L, lo, k, k_exp, congr) {

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

        if (lo > hi) {
            return nil
        }

        if (k == 1) {

            hi = max_p if (hi > max_p)
            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 (p.dec `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 t = lcm(L, p-1)
            m.is_coprime(t) || next

            __FUNC__(m*p, t, 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 3-Carmichael numbers in the range [1, 10^7] that are also strong pseudoprimes to base 2.

var k    = 3
var base = 2
var from = 1
var upto = 1e7

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

__END__
[15841, 29341, 52633, 252601, 314821, 1909001, 3581761, 4335241, 5049001, 5444489]