Starter code for LeetArxiv Semaev Index Calculus on secp elliptic curve experiment

SemaevIndexCalculusStarter.c

//Complete walkthrough: https://leetarxiv.substack.com/p/semaev-naive-index-calculus
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <stdbool.h>
#include <time.h>
#include <math.h>
//clear && gcc Semaev.c -lm -o m.o && ./m.o

struct Point_point_struct
{
	int x;
	int y;
	int infinity;
};
typedef struct Point_point_struct Point[1];

struct curve_struct
{
	int primeNumber;
	int pointOrder;
	int generatorX;
	int generatorY;
	int primitiveRoot;
	int b;
};

typedef struct curve_struct Curve[1];

void InitializeSmallCurve(Curve curve, int primeNumber, int pointOrder, int generatorX, int generatorY, int primitiveRoot)
{
	assert(primeNumber % 3 == 1);
	curve[0].primeNumber = primeNumber;
	curve[0].pointOrder = pointOrder;
	curve[0].generatorX = generatorX;
	curve[0].generatorY = generatorY;
	curve[0].primitiveRoot = primitiveRoot;
	curve[0].b = 7;
}

void PrintSmallCurve(Curve curve)
{
	printf("Prime number: %4d\n", curve->primeNumber);
	printf("Generator : (%4d, %4d)\n", curve->generatorX, curve->generatorY);
	printf("Point order: %4d\n", curve->pointOrder);
}

void InitializePoint(Point point)
{
	point[0].x = 0;
	point[0].y = 0;
	point[0].infinity = 0;
}

void CopyPoint(Point source, Point destination)
{
	destination->x = source->x;
	destination->y = source->y;
	destination->infinity = source->infinity;
}

void PrintPoint(Point point)
{
	printf("(%d %d)\n", point->x,point->y);
}

//Modular inverse of a mod m
int ModularInverse(int a, int m) 
{
	int t = 0, newT = 1;
	int r = m, newR = a;

	while(newR != 0)
	{
		int quotient = r / newR;
		int tempT = t;
		t = newT;
		newT = tempT - quotient * newT;
		int tempR = r;
		r = newR;
		newR = tempR - quotient * newR;
	}

	if(r > 1) return -1;  
	if(t < 0) t += m;

	return t;
}

bool PointIsOnCurve(Point point, int primeNumber)
{
	if(point->infinity) return true;
	int y2 = (point->y * point->y) % primeNumber;
	int x2 = (point->x * point->x) % primeNumber;
	int x3 = (x2 * point->x) % primeNumber;
	return y2 == (x3 + 7) % primeNumber;
}



bool AddCurvePoints(Point R, Point P, Point Q, int primeNumber)
{
	//Case 0: Handle Points at infinity
	if(P->infinity == 1){CopyPoint(Q,R);return true;}
	if(Q->infinity == 1){CopyPoint(P,R);return true;}
	
	//Case 1: Handle P->x == Q->x
	if(P->x == Q->x) 
	{
		//Case 1.1: X values similar but Y differ or 0
		if(P->y != Q->y || P->y == 0){R->infinity = 1;return true;}
		//Case 1.2: Point Doubling
		int num = P->x * P->x * 3;num = (num + primeNumber) % primeNumber;
		int den = P->y * 2;den = (den + primeNumber) % primeNumber;
		int denominatorInverse = ModularInverse(den, primeNumber);
		if(denominatorInverse == -1)
		{
			R->infinity = 1;
			return true;
		}
		//s = (3x^2)/(2y)
		int s = (num * denominatorInverse) % primeNumber;s = (s + primeNumber) % primeNumber;
		R->x = ((s*s) - 2 * P->x) % primeNumber;
		R->y = (s * (P->x - R->x) - P->y) % primeNumber;
		R->infinity = 0;
	}
	//Case 2: Handle P != Q (Point Addition)
	else
	{
		int num = Q->y - P->y;num = (num + primeNumber) % primeNumber;
		int den = Q->x - P->x;den = (den + primeNumber) % primeNumber;
		int denominatorInverse = ModularInverse(den, primeNumber);
		denominatorInverse = (denominatorInverse + primeNumber) % primeNumber;
		if(denominatorInverse == -1)
		{
			R->infinity = 1;
			return true;
		}
		int s = (num * denominatorInverse) % primeNumber;s = (s + primeNumber) % primeNumber;
		R->x = (s *s - P->x - Q->x) % primeNumber;
		//y3 = s*(x1 - x3) - y1
		R->y = (s *(P->x - R->x) - P->y) % primeNumber;
		R->infinity = 0;	
	}
	R->x = (R->x + primeNumber) % primeNumber;
	R->y = (R->y + primeNumber) % primeNumber;
	return true;
}

void LSB_DoubleAndAdd(Curve curve, Point result, int privateKey)
{
	//Create pointAtInfinity, temp0, temp1
	Point pointAtInfinity;
	InitializePoint(pointAtInfinity);
	pointAtInfinity->infinity = 1;  

	//Save generator to temp0
	Point temp0;
	InitializePoint(temp0);
	temp0->x = curve->generatorX;
	temp0->y = curve->generatorY;

	Point temp1;
	InitializePoint(temp1);
	assert(privateKey > 0);
	size_t binaryLength = floor(log2(privateKey)) + 1;
	for(size_t i = 0; i < binaryLength; ++i)
	{
		int bit = (privateKey >> i) & 1;
		if(bit)
		{
			//temp1 = pointAtInfinity + currentGenerator
			AddCurvePoints(temp1, pointAtInfinity, temp0, curve->primeNumber);
			CopyPoint(temp1, pointAtInfinity);
		}
		//Update temp0(currentGenerator) = 2*temp0(previousGenerator)
		AddCurvePoints(temp1, temp0, temp0, curve->primeNumber);
		CopyPoint(temp1, temp0);	
	}
	CopyPoint(pointAtInfinity, result);
}

void TestDoubleandAdd()
{
	int primeNumber = 20959;
	int pointOrder  = 20947;
	int generatorX  = 17263;
	int generatorY  = 13626;
	int primitiveRoot = 3;
	int b = 7;

	Curve smallCurve;
	InitializeSmallCurve(smallCurve,primeNumber,pointOrder,generatorX,generatorY,primitiveRoot);
	PrintSmallCurve(smallCurve);
	
	Point result;
	for(int privateKey = 1; privateKey < pointOrder; privateKey++)
	{
		LSB_DoubleAndAdd(smallCurve, result, privateKey);
		assert(PointIsOnCurve(result, smallCurve->primeNumber));
		printf("%4d: ",privateKey);PrintPoint(result);
	}
}

int main()
{
	TestDoubleandAdd();
	return 0;
}

StarterCodePython.py

#Complete walkthrough: https://leetarxiv.substack.com/p/semaev-naive-index-calculus
class Curve:
    def __init__(self, p, a, b, generator):
        self.p = p  # Prime modulus
        self.a = a  # Curve coefficient a
        self.b = b  # Curve coefficient b
        self.generator = Point(generator[0], generator[1], self)
    
    def is_on_curve(self, point):
        if point.infinity:
            return True
        x, y = point.x, point.y
        return (y * y - (x * x * x + self.a * x + self.b)) % self.p == 0

    def __repr__(self):
        return f"Curve over F_{self.p}: y^2 = x^3 + {self.a}x + {self.b}"


class Point:
    def __init__(self, x=None, y=None, curve=None, infinity=False):
        self.x = x
        self.y = y
        self.curve = curve
        self.infinity = infinity
    
    def __eq__(self, other):
        if self.infinity and other.infinity:
            return True
        if self.infinity or other.infinity:
            return False
        return self.x == other.x and self.y == other.y and self.curve == other.curve

    def __neg__(self):
        if self.infinity:
            return self
        return Point(self.x, (-self.y) % self.curve.p, self.curve)

    def __add__(self, Q):
        P = self
        if P.infinity:
            return Q
        if Q.infinity:
            return P

        if P.x == Q.x and (P.y != Q.y or P.y == 0):
            return Point(infinity=True, curve=self.curve)

        p = self.curve.p

        if P.x == Q.x:
            # Point doubling
            s_num = (3 * P.x * P.x + self.curve.a) % p
            s_den = modinv(2 * P.y, p)
        else:
            # Point addition
            s_num = (Q.y - P.y) % p
            s_den = modinv((Q.x - P.x) % p, p)

        s = (s_num * s_den) % p
        x_r = (s * s - P.x - Q.x) % p
        y_r = (s * (P.x - x_r) - P.y) % p
        return Point(x_r, y_r, self.curve)

    def __rmul__(self, k):
        result = Point(infinity=True, curve=self.curve)
        addend = self

        while k > 0:
            if k & 1:
                result = result + addend
            addend = addend + addend
            k >>= 1
        return result

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


def modinv(a, m):
    # Extended Euclidean Algorithm
    t, new_t = 0, 1
    r, new_r = m, a % m

    while new_r != 0:
        quotient = r // new_r
        t, new_t = new_t, t - quotient * new_t
        r, new_r = new_r, r - quotient * new_r

    if r != 1:
        raise ValueError(f"No inverse exists for {a} mod {m}")
    return t % m


# ======= Example Usage =======

if __name__ == "__main__":
    p = 20959
    a = 0
    b = 7
    Gx, Gy = 17263, 13626

    curve = Curve(p, a, b, generator=(Gx, Gy))
    G = curve.generator

    print("Elliptic Curve:", curve)
    print("Generator G:", G)
    print("Is G on the curve?", curve.is_on_curve(G))

    # Test scalar multiplication
    priv_key = 20942
    pub_key = priv_key * G

    print(f"{priv_key} * G = {pub_key}")
    print("Is pub_key on curve?", curve.is_on_curve(pub_key))