OEIS strong elliptic Carmichael numbers.sagews

def exponent(p,power,E):
    #Assumes p does not divide |E(F_p)|.
    #For the curve E:y^2 = x^3 + 80, this gives the correct result because for all primes p, p does not divide |E(F_p)|.
    Ep = E.reduction(p)
    exp = Ep.abelian_group().generator_orders()[0]
    return exp*p^(power-1) 
def is_strong_elliptic_Carmichael_number(E,N):
    if gcd(E.discriminant(),N) > 1:
        return false
    factorization = list(factor(N))
    if len(factorization) <= 1: return false
    cond = E.conductor()
    aN = 1
    prime_coefficients = {}
    #compute ap, aN
    for prime, power in factorization:
        listy = [1, prime+1-E.Np(prime)]
        char = 1
        if mod(cond,prime) == 0:
            char = 0
        for j in range(2,power+1):
            listy.append(listy[1]*listy[j-1]-char*prime*listy[j-2])
        aN = aN*listy[-1]
        prime_coefficients[prime] = listy[1]
    t = N+1-aN
    while mod(t,2) != 1:
        t/=2
    #check exponent of group E(Z/p^ord_p(N)Z) divides t.
    for prime, power in factorization:
        exp = exponent(prime,power,E)
        if mod(t, exp) != 0:
            return false
    return true

A,B = 0,80
E = EllipticCurve([A,B])
count = 1
N = 1
while true:
    if is_strong_elliptic_Carmichael_number(E,N):
        print count,N
        count+=1
    N+=1

1 7783 2 327781 3 416521 4 991681 5 3401341 6 4584187 7 5234089 8 5325883 9 991681 5 3401341 6 4584187 7 5234089 8 5325883 9