experimental-projects/pseudoprimes/generators/imprimitive_carmichael_coin_change.sf

#!/usr/bin/ruby

# Generate imprimitive Carmichael numbers with n prime factors.

func is_imprimitive(n) {
    n.factor.gcd_by { .dec }**2 > lambda(n)
}

var denominations = []

func change(m, pos=0, solution=[]) {

    if (solution.len >= 2) {
        if (solution[-1] / solution[0] > 40) {
            #say '# here'
            return nil
        }
    }

    if (solution.len == m) {

        with (solution.prod) {|C|
            if (C.is_carmichael) {
                say C
                say "# Imprimitive: #{C}" if is_imprimitive(C)
            }
        }

        return nil
    }

    if (solution.len > m) {
        return nil
    }

    if (pos > denominations.end) {
        return nil;
    }

    change(m, pos + 1, solution);
    change(m, pos + 1, [solution..., denominations[pos]]);
}

func generate_imprimitive(p, m) {

    for z in (2..100) {

        var arr = []

        for k in (1 .. 2.sqrt**z) {
            k.is_smooth(5) || next
            #next if (2*5 > 40)
            var r = (2*k*p + 1)
            if (r.is_prime && (r.dec.gpf == p)) {
                arr << r
            }
        }

        var t = binomial(arr.len, m)

        t >= m  || next
        t < 1e6 || break

        arr = arr.last(19448)

        var count = 0
        say "# Combinations: #{t}"

        denominations = arr
        change(m)

        #~ arr.combinations(m, {|*a|
            #~ with (a.prod) { |C|
                #~ break if (C > 325533792014488126487416882038879701391121)
                #~ #break if (C.gpf / C.lpf > 40)
                #~ if (C.is_carmichael) {
                    #~ say C
                    #~ say "# Imprimitive with p = #{p}: #{C}" if is_imprimitive(C)
                #~ }
            #~ }
            #~ #break if (++count > 1e5)
        #~ })
    }
}

for p in (primes(73..457)) {
     say "# p = #{p}"
     generate_imprimitive(p,10)
}