sidef-scripts/Math/strong_fermat_pseudoprimes_in_range.sf

#!/usr/bin/ruby

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

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

#`(
# PARI/GP program (slow):

strong_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, p, j, k_exp, congr) = my(list=List()); forprime(q=p, sqrtnint(B\m, j), if(base%q != 0, my(tv=valuation(q-1, 2)); if(tv > k_exp && Mod(base, q)^(((q-1)>>tv)<<k_exp) == congr, my(v=m*q, t=q, r=nextprime(q+1)); while(v <= B, my(L=lcm(l, znorder(Mod(base, t)))); if(gcd(L, v) == 1, if(j==1, if(v>=A && if(k==1, !isprime(v), 1) && (v-1)%L == 0, listput(list, v)), if(v*r <= B, list=concat(list, f(v, L, r, j-1, k_exp, congr)))), break); v *= q; t *= q)))); list); my(r=f(1, 1, 2, k, 0, 1)); for(v=0, logint(B, 2), r=concat(r, f(1, 1, 2, k, v, -1))); vecsort(Vec(r));

# PARI/GP program (fast):

strong_check(p, base, e, r) = my(tv=valuation(p-1, 2)); tv > e && Mod(base, p)^((p-1)>>(tv-e)) == r;
strong_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k, e, r) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, forstep(p=lift(1/Mod(m, l)), hi, l, if(isprimepower(p) && gcd(m*base, p) == 1 && strong_check(p, base, e, r), my(n=m*p); if(n >= A && (n-1) % znorder(Mod(base, p)) == 0, listput(list, n)))), forprime(p=lo, hi, base%p == 0 && next; strong_check(p, base, e, r) || next; my(z=znorder(Mod(base, p))); gcd(m,z) == 1 || next; my(q=p, v=m*p); while(v <= B, list=concat(list, f(v, lcm(l, z), p+1, k-1, e, r)); q *= p; Mod(base, q)^z == 1 || break; v *= p))); list); my(res=f(1, 1, 2, k, 0, 1)); for(v=0, logint(B, 2), res=concat(res, f(1, 1, 2, k, v, -1))); vecsort(Set(res));

)

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

    A = max(k.pn_primorial, A)

    var seen = Hash()

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

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

        lo > hi && return nil

        if (j == 1) {

            if (L == 1) {   # optimization
                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

                    for (var v = (m == 1 ? p*p : m*p); v <= B; v *= p) {
                        v >= A                || next
                        v.is_strong_psp(base) || break
                        callback(v) if !(seen{v} := 0 ++)
                    }
                })
                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_power   || next
                p.is_coprime(m)    || next
                p.is_coprime(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 v = m*p
                v >= A || next

                if (znorder(base, p) `divides` v.dec) {
                    callback(v) if !(seen{v} := 0 ++)
                }
            }

            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

            var q = p
            var v = m*p

            while (v <= B) {
                __FUNC__(v, lcm(L, z), p+1, j-1, k_exp, congr)
                q *= p
                powmod(base, z, q) == 1 || break
                v *= p
            }
        })
    }

    # 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 strong Fermat pseudoprimes to base 3 in range [1, 10^4]

var from = 1
var upto = 1e4
var base = 3
var arr = []

for k in (1..100) {
    break if (k.pn_primorial > upto)
    strong_fermat_pseudoprimes_in_range(from, upto, k, base, { arr << _ })
}

say arr.sort

# Run some tests

if (false) {    # true to run some tests
    for k in (1..5) {

        say "Testing k = #{k}"

        var lo = k.pn_primorial
        var hi = lo*1000

        for base in (2..100) {

            var this = k.strong_fermat_psp(base, lo, hi)
            #var this = k.omega_primes(lo,hi).grep{.is_strong_psp(base) && .is_composite }
            var that = gather {
                strong_fermat_pseudoprimes_in_range(lo, hi, k, base, func (n) { take(n) })
            }.sort

            this == that ||
                die "Error for k = #{k} and base = #{base} with hi = #{hi}\n#{this} != #{that}";
        }
    }
}

__END__
[121, 703, 1891, 3281, 8401, 8911]