sidef-scripts/Math/smallest_carmichael_divisible_by_n_faster.sf

#!/usr/bin/ruby

# Method for finding the smallest Carmichael number divisible by n.

# See also:
#   https://oeis.org/A135721
#   https://oeis.org/A253595

func carmichael_numbers_from_multiple(A, B, m, L, lo, k, callback) {

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

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

    return nil if (lo > hi)

    if (k == 1) {

        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
            p.divides(m) && next
            with (m*p) {|n|
                if (p.dec `divides` n.dec) {
                    callback(n)
                }
            }
        }

        return nil
    }

    each_prime(lo, hi, {|p|

        p.divides(m) && next

        var t = lcm(L, p-1)
        m.is_coprime(t) || next

        __FUNC__(A, B, m*p, t, p+1, k-1, callback)
    })
}

func carmichael_divisible_by(m) {

    m >= 1 || return nil
    m.is_even && return nil
    gcd(m, phi(m)) == 1 || return nil

    var a = max(561, m)
    var b = 2*a
    var L = m.factor.lcm{.dec}

    var found = []

    loop {

        for k in ((m.is_prime ? 2 : 1)..1000) {

            var P = k.by {|p|
                p.is_odd && p.is_prime && !p.divides(m) && !p.divides(L)
            }

            break if (P.prod*m > b)

            var callback = {|n|
                found << n
                b = min(b, n)
            }

            carmichael_numbers_from_multiple(a, b, m, L, P[0], k, callback)
        }

        a = b+1
        b = 2*a

        break if found
    }

    found.min
}

assert_eq(carmichael_divisible_by(3), 561)
assert_eq(carmichael_divisible_by(3*5), 62745)
assert_eq(carmichael_divisible_by(7*19), 1729)
assert_eq(carmichael_divisible_by(47*89), 62745)
assert_eq(carmichael_divisible_by(5*47*89), 62745)
assert_eq(carmichael_divisible_by(3*47*89), 62745)
assert_eq(carmichael_divisible_by(3*89), 62745)

say carmichael_divisible_by.map(primes(3..50))
say 40.of(carmichael_divisible_by).grep

__END__
[561, 1105, 1729, 561, 1105, 561, 1729, 6601, 2465, 2821, 29341, 6601, 334153, 62745]
[561, 561, 1105, 1729, 561, 1105, 62745, 561, 1729, 6601, 2465, 2821, 561, 825265, 29341]