Endomorphism structure of y^2=x^3+3 % 199

README.md

Uses Plotly and igraph to render an interactive 3d graph of point operations on a j-invariant 0 elliptic curve defined over p=199 via y^2+x^3+3.

endomorphine-structure.py

import sys
import math
import plotly.graph_objects as go
import igraph as ig
from collections import namedtuple
from dataclasses import dataclass
import z3

Node = namedtuple('Node',['name','group'])
GraphEdgesT = set[tuple[str,str]]
GraphNodesT = dict[str,Node]

@dataclass
class Node:
    name:str
    group:str

@dataclass
class GraphT:
    nodes: GraphNodesT
    edges: GraphEdgesT
    def add(self,node_name:str,color:str,targets:list[tuple[str,str]]=[]):
        if node_name not in self.nodes:
            self.nodes[node_name] = Node(node_name, color)
        for t in targets:
            self.edges.add((node_name,t))
    def __add__(self, other:'GraphT') -> 'GraphT':
        return graph_combine(self, other)

def graph_combine(*args:GraphT):
    nodes, edges = dict(),set()
    for a in args:
        nodes.update(a.nodes)
        edges = edges.union(a.edges)
    return GraphT(nodes,edges)

def graph_linearize(args:GraphT):
    nodes = args.nodes
    name_to_index = {name: i for i, name in enumerate(nodes.keys())}
    # Convert edges to use numeric indices
    edges = [(name_to_index[source], name_to_index[target])
            for source,target in list(sorted(args.edges))]
    labels = [node.name for node in nodes.values()]
    groups = [node.group for node in nodes.values()]
    return (nodes, labels, groups), edges

def flatten(xss):
    return [x for xs in xss for x in xs]

def squares(x, p):
    while True:
        x = (x * x) % p
        yield x

###############################################################################

class Point(namedtuple('_Point', ['x','y','a','b','p','is_infinity'])):
    """Class representing a point on an elliptic curve."""
    @classmethod
    def infinity(cls, a, b, p):
        """Return the point at infinity."""
        return cls(0, 0, a, b, p, True)

    @property
    def coords(self):
        return (self.x, self.y)

    @property
    def curve(self):
        return (self.a, self.b, self.p)

    def is_on_curve(self):
        """Check if the point lies on the curve."""
        if self.is_infinity:
            return True
        left = (self.y * self.y) % self.p
        right = (self.x * self.x * self.x + self.a * self.x + self.b) % self.p
        return left == right

    def __eq__(self, other:'Point'):
        if self.curve != other.curve:
            return False
        if self.is_infinity:
            return other.is_infinity
        if other.is_infinity:
            return False
        return self.x == other.x and self.y == other.y

    def __str__(self):
        if self.is_infinity:
            return "Point(infinity)"
        return f"Point({self.x}, {self.y})"

    def __add__(self, other:'Point') -> 'Point':
        """Add two points using the elliptic curve group law."""
        assert self.curve == other.curve
        if self.is_infinity:
            return other
        if other.is_infinity:
            return self

        # Point doubling
        if self.x == other.x:
            if (self.y + other.y) % self.p == 0:
                return Point.infinity(self.a, self.b, self.p)
            else:
                # Compute the slope of the tangent line
                lambda_val = ((3 * self.x * self.x + self.a) * pow(2 * self.y, -1, self.p)) % self.p
        else:
            # Point addition
            lambda_val = ((other.y - self.y) * pow(other.x - self.x, -1, self.p)) % self.p

        x3 = (lambda_val * lambda_val - self.x - other.x) % self.p
        y3 = (lambda_val * (self.x - x3) - self.y) % self.p
        return Point(x3, y3, self.a, self.b, self.p, False)

    def scalar_mul(self, k):
        """Multiply point by scalar k using double-and-add algorithm."""
        result = Point.infinity(self.a, self.b, self.p)
        addend = self
        while k > 0:
            if k & 1:
                result = result + addend
            addend = addend + addend
            k >>= 1
        return result

###############################################################################

def draw_graph(output,graph:GraphT,modifiers_by_color:dict[tuple[int,int]]):
    (nodes, labels, groups), edges = graph_linearize(graph)

    N = len(nodes)
    G = ig.Graph(edges, directed=False)

    # Set vertex names as attributes
    G.vs['name'] = labels
    G.vs['group'] = groups

    layt = G.layout('kk', dim=3, kkconst=N*5)
    for i in range(len(layt)):
        color = groups[i]
        if color in modifiers_by_color:
            (zoffset, multiplier) = modifiers_by_color[color]
            layt[i][2] = layt[i][2] + zoffset
            layt[i][0] *= multiplier
            layt[i][1] *= multiplier
            layt[i][2] *= multiplier

    Xn = [layt[k][0] for k in range(N)]
    Yn = [layt[k][1] for k in range(N)]
    Zn = [layt[k][2] for k in range(N)]

    Xe = flatten([[layt[e[0]][0], layt[e[1]][0], None] for e in edges])
    Ye = flatten([[layt[e[0]][1], layt[e[1]][1], None] for e in edges])
    Ze = flatten([[layt[e[0]][2], layt[e[1]][2], None] for e in edges])

    trace1 = go.Scatter3d(x=Xe,
                y=Ye,
                z=Ze,
                mode='lines',
                line=dict(color='rgb(125,125,125)', width=1),
                hoverinfo='none')

    trace2 = go.Scatter3d(x=Xn,
                y=Yn,
                z=Zn,
                mode='markers',
                name='actors',
                marker=dict(symbol='circle',
                            size=6,
                            color=groups,
                            colorscale='Viridis',
                            line=dict(color='rgb(50,50,50)',width=0.5)),
                text=labels,
                hoverinfo='text')

    axis = dict(showbackground=False,
            showline=False,
            zeroline=False,
            showgrid=False,
            showticklabels=False,
            title='')

    layout = go.Layout(
            title="Endomorphism structure of y^2=x^3+3 % 199",
            width=1000,
            height=1000,
            showlegend=False,
            scene=dict(xaxis=dict(axis), yaxis=dict(axis), zaxis=dict(axis)),
        margin=dict(t=100),
        hovermode='closest')

    fig = go.Figure(data=[trace1, trace2], layout=layout)
    fig.write_html(output)

###############################################################################

def squarepairs(x, p):
    prev = x
    for nxt in squares(x, p):
        yield (str(prev),str(nxt))
        prev = nxt

def squaregraph(p, color:str):
    g = GraphT(dict(), set())
    for i in range(1,p):
        g.nodes[str(i)] = Node(str(i), color)
        for k in squarepairs(i, p):
            if k in g.edges:
                break
            g.edges.add(k)
    return g

def mulgraph(p,n):
    g = GraphT(dict(), set())
    for a in range(1,n):
        for b in range(1,n):
            if a != b:
                g.add(f'{a}+{b}', 'orange', [str(a), str(b), str((a+b)%p)])
    return g

################################################################################

def find_double_chain(p:Point):
    g = GraphT(dict(), set())
    current = p
    while True:
        next:Point = current + current
        name = f'E{p.coords}*2'
        if name in g.nodes:
            break
        x,y = current.coords
        g.add(str(current.coords), 'blue', [f'x={x}', f'y={y}'])
        g.add(f'x={x}', 'orange', [f'y={y}'])
        g.add(f'y={y}', 'yellow', [f'x={x}'])
        g.nodes[name] = Node(name, 'magenta')
        g.edges.add((str(p.coords),name))
        g.edges.add((name,str(next.coords)))
        p = next
    return g

def find_triple_chain(p:Point):
    g = GraphT(dict(), set())
    current = p
    while True:
        next:Point = current + current + current
        name = f'E{p.coords}*3'
        if name in g.nodes:
            break
        x,y = current.coords
        g.add(str(current.coords), 'blue', [f'x={x}', f'y={y}'])
        g.add(f'x={x}', 'orange', [f'y={y}'])
        g.add(f'y={y}', 'yellow', [f'x={x}'])
        g.nodes[name] = Node(name, 'magenta')
        g.edges.add((str(p.coords),name))
        g.edges.add((name,str(next.coords)))
        p = next
    return g

def find_quad_chain(p:Point):
    g = GraphT(dict(), set())
    current = p
    while True:
        next:Point = current + current + current + current
        name = f'E{p.coords}*4'
        if name in g.nodes:
            break
        x,y = current.coords
        g.add(str(current.coords), 'blue', [f'x={x}', f'y={y}'])
        g.add(f'x={x}', 'orange', [f'y={y}'])
        g.add(f'y={y}', 'yellow', [f'x={x}'])
        g.nodes[name] = Node(name, 'deeppink')
        g.edges.add((str(p.coords),name))
        g.edges.add((name,str(next.coords)))
        p = next
    return g

from typing import Callable, Optional
PointFilterT = Callable[[int,int],bool]

def curvepoints(p, a, b,filter:Optional[PointFilterT]=None):
    g = GraphT(dict(), set())
    #g.add(str((0,0)), 'red')
    g.add(str((0,0)),"red", ['infinity'])
    for x in range(p):
        for y in range(p):
            point = Point(x, y, a, b, p, False)
            if not point.is_on_curve():
                continue
            if filter and not filter(x,y):
                continue
            name = str(point.coords)
            #g.add(name, 'blue')
            #g.add(name, 'blue', [f'{x}', f'{y}'])
            g.add(name, 'blue', [f'x={x}', f'y={y}'])
            g.add(f'x={x}', 'orange', [f'y={y}'])
            g.add(f'y={y}', 'yellow', [f'x={x}'])
            if point.is_infinity:
                g.nodes['infinity'] = Node('infinity','red')
                g.edges.add(('infinity',name))
            g += find_double_chain(point)
            #g += find_oct_chain(point)
            #g += find_triple_chain(point)
            #g += find_quad_chain(point)
            #g += find_addition_chain(point)
    return g

###############################################################################

def find_endomorphisms(p,a,b,filter:Optional[PointFilterT]=None):
    g = GraphT(dict(),set())
    all_k = set()
    for x in range(p):
        for y in range(p):
            point = Point(x,y,a,b, p, x==0 or y==0)
            if filter and not filter(x,y):
                continue
            if point.is_infinity or not point.is_on_curve():
                continue
            #g.add(str(point.coords), 'purple')
            for phi_x in range(2,p):
                looped = set()
                k =f'Φx({phi_x})@'
                prev = None
                while True:
                    x2 = (point.x*phi_x)%p
                    p2 = Point(x2,point.y,point.a, point.b, point.p,x2 == 0 or point.is_infinity)
                    point_coords = str(point.coords)
                    combined = ''.join([k, point_coords])
                    if combined in looped:
                        break
                    if p2.coords == point.coords:
                        looped.add(combined)
                        continue
                    if p2.is_infinity or not p2.is_on_curve():
                        looped.add(combined)
                        continue

                    looped.add(combined)
                    g.add(combined, 'grey')
                    if prev is not None:
                        prev_combined = ''.join([k, prev])
                        g.edges.add((prev_combined, combined))
                        prev = p2.coords
                        print('Loop')

                    next_coords = str(p2.coords)
                    #g.edges.add((point_coords, combined))
                    g.edges.add((f'x={point.x}', combined))
                    g.edges.add((f'y={point.y}', combined))
                    g.edges.add((combined, next_coords))

            for phi_y in range(2,p):
                looped = set()
                k =f'Φy({phi_y})@'
                prev = None
                while True:
                    y2 = (point.y*phi_y)%p
                    p2 = Point(point.x,y2,point.a, point.b, point.p,x2 == 0 or point.is_infinity)
                    point_coords = str(point.coords)
                    combined = ''.join([k, point_coords])
                    if combined in looped:
                        break
                    if p2.coords == point.coords:
                        looped.add(combined)
                        continue
                    if p2.is_infinity or not p2.is_on_curve():
                        looped.add(combined)
                        continue

                    looped.add(combined)
                    g.add(combined, 'pink')
                    if prev is not None:
                        prev_combined = ''.join([k, prev])
                        g.edges.add((prev_combined, combined))
                        prev = p2.coords
                        print('Loop')

                    next_coords = str(p2.coords)
                    #g.edges.add((point_coords, combined))
                    g.edges.add((f'x={point.x}', combined))
                    g.edges.add((f'y={point.y}', combined))
                    g.edges.add((combined, next_coords))

            #"""
            for phi_y in range(2,p):
                looped = set()
                k =f'Φ-x,y({phi_y})@'
                prev = None
                while True:
                    #x2 = point.x
                    #((p+1)/4)
                    x2 = (-point.x)%p
                    y2 = (point.y*phi_y)%p
                    #y2 = (-point.y)%p
                    #x2 = -(point.x)%p
                    #y2 = (point.y*phi_y)%p
                    #x2 = (point.x*point.x)%p
                    #y2 = (point.y*phi_y)%p
                    p2 = Point(x2,y2,point.a, point.b, point.p,x2 == 0 or point.is_infinity)
                    point_coords = str(point.coords)
                    combined = ''.join([k, point_coords])
                    if combined in looped:
                        break
                    if p2.coords == point.coords:
                        looped.add(combined)
                        continue
                    if p2.is_infinity or not p2.is_on_curve():
                        looped.add(combined)
                        continue

                    looped.add(combined)
                    g.add(combined, 'black')
                    if prev is not None:
                        prev_combined = ''.join([k, prev])
                        g.edges.add((prev_combined, combined))
                        prev = p2.coords
                        print('Loop')

                    next_coords = str(p2.coords)
                    #g.edges.add((point_coords, combined))
                    g.edges.add((f'x={point.x}', combined))
                    g.edges.add((f'y={point.y}', combined))
                    g.edges.add((combined, next_coords))
                    print("\t", point, f'Φ = (-x,{phi_y}*y) =', p2)
            #"""

            """
            for phi_x in range(2,p):
                phi_y = phi_x
                #for phi_y in range(2,p):
                looped = set()
                k =f'Φx(^{phi_x})y(^{phi_y})@'
                prev = None
                while True:
                    x2 = pow(point.x,phi_x,p)
                    y2 = pow(point.y,phi_y,p)
                    p2 = Point(x2,y2,point.a, point.b, point.p,x2 == 0 or point.is_infinity)
                    point_coords = str(point.coords)
                    combined = ''.join([k, point_coords])
                    if combined in looped:
                        break
                    if p2.coords == point.coords:
                        looped.add(combined)
                        continue
                    if p2.is_infinity or not p2.is_on_curve():
                        looped.add(combined)
                        continue

                    looped.add(combined)
                    g.add(combined, 'lightseagreen')
                    if prev is not None:
                        prev_combined = ''.join([k, prev])
                        g.edges.add((prev_combined, combined))
                        prev = p2.coords
                        print('Loop')

                    all_k.add(k)
                    next_coords = str(p2.coords)
                    g.edges.add((f'x={point.x}', combined))
                    g.edges.add((f'y={point.y}', combined))
                    g.edges.add((combined, next_coords))
            #"""
    print(all_k)

    return g

###############################################################################

def main(p,g):
    p = int(p)
    g = int(g)
    graph = GraphT(dict(),set())
    a = 0
    q = 211
    b = 3**1
    filter = lambda x,y: Point(x,y,a,b,p,False).is_on_curve()
    #graph += squaregraph(p,'aquamarine')
    for i in [1,5]: # range(6):
        graph += curvepoints(p, a, pow(b,i,p), filter)
        graph += find_endomorphisms(p,a,pow(b,i,p), filter)

    g = Point(182,142,a,b,p,False)
    start = g
    i = 1
    graph.add('infinity', 'black')
    graph.add('0', 'black', ['infinity'])
    while True:
        print("Point", i, start)
        k = str(start.coords)
        if k not in graph.nodes:
            graph.nodes[k] = Node(k, 'blue')
        graph.add(str(i),"darkorchid", [k])
        #next_i = (i+1)%q
        #next_i = pow(2,(p+1)//8,q)
        graph.edges.add((str(i), str((i+i)%q)))
        #graph.edges.add((str(i), str(pow(i,2,q))))
        graph.edges.add((str(i), str((-i)%q)))
        #if i > 1:
        #    graph.add(str((i-1)),"darkorchid", [str(i)])
        start += g
        if start.is_infinity or not start.is_on_curve():
            break
        i += 1
        if g == start:
            break

    draw_graph("endomorphine-structure.html", graph, {
        #'purple': (0, 10),
    })
    print("p%12",p%12)

###############################################################################

if __name__ == "__main__":
    sys.exit(main(*sys.argv[1:]))