carmichael_numbers_validity.py

from matplotlib import pyplot as plt

# http://oeis.org/wiki/Carmichael_numbers
CARMICHAEL_NUMBERS = [
    561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657,
    52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461,
    252601, 278545, 294409, 314821, 334153
]

def prime_factors(n):
    """Get the prime factors of n."""
    i = 2
    factors = []
    while i * i <= n:
        if n % i:
            i += 1
        else:
            n //= i
            factors.append(i)
    if n > 1:
        factors.append(n)
    return factors

if __name__ == '__main__':
    percentages = []

    for carmichael_number in CARMICHAEL_NUMBERS:
        progressions_to_test = carmichael_number - 1
        factors = prime_factors(carmichael_number)

        valid_progressions = []

        for i in range(1, progressions_to_test + 1):
            if [i % factors[j] != 0 for j in range(len(factors))] == [True] * len(factors):
                valid_progressions.append(i)

        print('For {}, {}% of the first 1000 progressions verify the property'.format(
            carmichael_number,
            len(valid_progressions) / progressions_to_test * 100))

        percentages.append(len(valid_progressions) / progressions_to_test)

    plt.plot(percentages)
    plt.show()