for E in cremona_optimal_curves([389..30000]):
    r = E.rank(use_database=True)
    if r >= 2:
        N = E.conductor()
        for F in cremona_optimal_curves([N]):
            if F != E:
                for p in prime_divisors(gcd(E.modular_degree(), F.modular_degree())):
                    if p > 2:
                        if all( (E.ap(ell)-F.ap(ell))%p == 0 for ell in primes(2*N) if (N*p) % ell ):
                            if E.galois_representation().is_surjective(p):
                                print E.cremona_label(), p, F.cremona_label()
                                sys.stdout.flush()