Verification script for https://twitter.com/pbarreto/status/869103226276134912

Curve13318.sage

#!/usr/bin/env sage
# *-* encoding: utf-8 *-*

"""
Find a value `b` for the curve `E : y^2 = x^3 - 3x + b` s.t.
the curve has a prime order and its twist also has a prime order.
"""

from itertools import chain, count

F = GF(2^255 - 19)

def plusmin(it):
    for x in it:
        yield x
        yield -x

for i,b in enumerate(plusmin(count())):
    percents = int(100.0 * i / float(2*13318))
    print('[% 3d%%] b = % 6d' % (percents, b))

    try:
        E = EllipticCurve(F, [-3, b])
    except ArithmeticError:
        # Curve probably has a singularity
        continue
    n = E.order()

    # Only accept curves of prime order
    if not is_prime(n):
        continue

    # Only accept if the twist is also of prime order
    twist = E.quadratic_twist()
    if not is_prime(twist.order()):
        continue

    bits = int(log(float(E.order()))/log(2.0))
    print('')
    print('Found a good curve E : y^2 = x^3 - 3x + {}'.format(b))

    # Start selecting random points for the rest of the search to find G
    for x in plusmin(range(0, 256)):
        for y in range(256):
            try:
                G = E(x, y)
            except TypeError: # (x, y) is not a valid point
                continue
            print('  - G = {}'.format(G))
    break

Curve_aZero.sage

#!/usr/bin/env sage
# *-* encoding: utf-8 *-*

"""
Find a value `b` for the curve `E : y^2 = x^3 - 3x + b` s.t.
the curve has a prime order and its twist also has a prime order.
"""

from itertools import chain, count

F = GF(2^255 - 19)

def plusmin(it):
    for x in it:
        yield x
        yield -x

for i,b in enumerate(plusmin(count())):
    percents = int(100.0 * i / float(2*13318))
    print('[% 3s%%] b = % 6d' % ("??", b))

    try:
        E = EllipticCurve(F, [0, b])
    except ArithmeticError:
        # Curve probably has a singularity
        continue
    n = E.order()

    # Only accept curves of prime order
    if not is_prime(n):
        continue

    # Only accept if the twist is also of prime order
    twist = E.quadratic_twist()
    if not is_prime(twist.order()):
        continue

    bits = int(log(float(E.order()))/log(2.0))
    print('')
    print('Found a good curve E : y^2 = x^3 - 3x + {}'.format(b))

    # Start selecting random points for the rest of the search to find G
    for x in plusmin(range(0, 256)):
        for y in range(256):
            try:
                G = E(x, y)
            except TypeError: # (x, y) is not a valid point
                continue
            print('  - G = {}'.format(G))
    break

Runtime is around 24 hours.