sidef-scripts/Math/smallest_lucas-carmichael_divisible_by_n_faster.sf

#!/usr/bin/ruby

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

# See also:
#   https://oeis.org/A253597
#   https://oeis.org/A253598

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

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

    return nil if (lo > hi)

    if (k == 1) {

        lo = max(lo, idiv_ceil(A, m))
        lo > hi && return nil

        var t = mulmod(m.invmod(L), -1, 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.inc `divides` n.inc) {
                    callback(n)
                }
            }
        }

        return nil
    }

    each_prime(lo, hi, {|p|

        p.divides(m) && next
        m.is_coprime(p+1) || next

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

func lucas_carmichael_divisible_by(m) {

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

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

    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)
            }

            lucas_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(lucas_carmichael_divisible_by(3), 399)
assert_eq(lucas_carmichael_divisible_by(3*7), 399)
assert_eq(lucas_carmichael_divisible_by(7*19), 399)

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

__END__
[399, 935, 399, 935, 2015, 935, 399, 4991, 51359, 2015, 1584599, 20705, 5719, 18095]
[399, 399, 935, 399, 935, 2015, 935, 399, 399, 4991, 51359, 2015, 8855, 1584599, 9486399]