$ a="486662"; b="1"; p="57896044618658097711785492504343953926634992332820282019728792003956564819949"
$ x="9"; y="43114425171068552920764898935933967039370386198203806730763910166200978582548"
$ ./ec.out "$a" "$b" "$p" "$x" "$y"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    25401056036235045220436215630227531782477232918721757282083167532085
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (2152181765955508802144811366391548358206596269149408838626058763387933371482, 11578758651574146923540708489424215527254544219548915855368449621649032022464)

$ ./ec.out "$a" "$b" "$p" "$x" "$y"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    14923829837058818803507170043209691418601607100463140519943405743260
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (10616717952671290844210554362351196131260739681824015303968067548318407877443, 44474405861832920486283648404346377112839096552743000067399554821542048806382)

$ x="10616717952671290844210554362351196131260739681824015303968067548318407877443"; y="44474405861832920486283648404346377112839096552743000067399554821542048806382"
$ ./ec.out "$a" "$b" "$p" "$x" "$y" "25401056036235045220436215630227531782477232918721757282083167532085"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    25401056036235045220436215630227531782477232918721757282083167532085
      * P(x, y) = (10616717952671290844210554362351196131260739681824015303968067548318407877443, 44474405861832920486283648404346377112839096552743000067399554821542048806382)

      = Q(x, y) = (5796465362698668629448044005371479448781495853974848652041497206752313461872, 30198337259834952114045673142150509214022566212693868468740722556827195886413)

$ x="2152181765955508802144811366391548358206596269149408838626058763387933371482"; y="11578758651574146923540708489424215527254544219548915855368449621649032022464"
$ ./ec.out "$a" "$b" "$p" "$x" "$y" "14923829837058818803507170043209691418601607100463140519943405743260"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    14923829837058818803507170043209691418601607100463140519943405743260
      * P(x, y) = (2152181765955508802144811366391548358206596269149408838626058763387933371482, 11578758651574146923540708489424215527254544219548915855368449621649032022464)

      = Q(x, y) = (5796465362698668629448044005371479448781495853974848652041497206752313461872, 30198337259834952114045673142150509214022566212693868468740722556827195886413)

$ 

#include "bn.c"

struct ecurve {
	bnum *a, *b, *p, *x, *y;
};

#define ecc struct ecurve

struct ectemp {
	bnum *i, *s, *xr, *yr;
	bnum *t, *u, *v;
	bnum *w, *h, *g;
};

#define ect struct ectemp

ecc *ecinit(bnum *a, bnum *b, bnum *p, bnum *x, bnum *y)
{
	ecc *r = malloc(1 * sizeof(ecc));
	r->a = a; r->b = b; r->p = p;
	r->x = x; r->y = y;
	return r;
}

ecc *ecdup(ecc *e)
{
	ecc *r = malloc(1 * sizeof(ecc));
	r->a = bndup(e->a); r->b = bndup(e->b); r->p = bndup(e->p);
	r->x = bndup(e->x); r->y = bndup(e->y);
	return r;
}

void ecfree(ecc *e)
{
	bnfree(e->a); bnfree(e->b); bnfree(e->p);
	bnfree(e->x); bnfree(e->y);
	free(e);
}

void ecout(int d, char *s, ecc *e, char *t)
{
	char *a = bnstr(e->a), *b = bnstr(e->b), *p = bnstr(e->p);
	char *x = bnstr(e->x), *y = bnstr(e->y);
	char as[2], bs[2];
	as[0] = '+'; as[1] = '\0';
	bs[0] = '\0'; bs[1] = '\0';
	if ((e->a)->sign == 1) { as[0] = '-'; }
	if ((e->b)->sign == 1) { bs[0] = '-'; }
	if (d == 1) { printf("  %s%s*(y^2) = x^3 %s %s*(x^2) + x (mod %s)\n", bs, b, as, a, p); }
	printf("%s", s);
	printf("(x, y) = (%s, %s)", x, y);
	printf("%s", t);
	free(a); free(b); free(p);
	free(x); free(y);
}

ect *etinit(ecc *e)
{
	ect *t = malloc(1 * sizeof(ect));
	int ss = max(1, (e->b)->size);
	ss = max(((e->a)->size * 2) + 2, ((e->p)->size * 2) + 2);
	ss = max(((e->x)->size * 2) + 2, ((e->y)->size * 2) + 2);
	t->i = bninit(ss); t->s = bninit(ss); t->xr = bninit(ss); t->yr = bninit(ss);
	t->t = bninit(ss); t->u = bninit(ss); t->v = bninit(ss);
	int tt = ((ss * 2) + 2);
	t->w = bninit(tt); t->h = bninit(tt); t->g = bninit(tt + 4);
	return t;
}

void etfree(ect *t)
{
	bnfree(t->i); bnfree(t->s); bnfree(t->xr); bnfree(t->yr);
	bnfree(t->t); bnfree(t->u); bnfree(t->v);
	bnfree(t->w); bnfree(t->h); bnfree(t->g);
	free(t);
}

// modular multiplicative inverse

void egcd(bnum *a, bnum *b, bnum *g)
{
	int size = ((a->size + b->size) * 3);
	// s = 0; news = 1
	bnum *s = bninit(size);
	bnum *news = bninit(size); news->nums[0] = 1;
	// r = b; newr = a
	bnum *r = bninit(size); bncopy(b, r);
	bnum *newr = bninit(size); bncopy(a, newr);
	// init some temp vars
	bnum *prev = bninit(size), *quot = bninit(size), *temp = bninit(size);
	while ((r->leng > 1) || (r->nums[0] > 0))
	{
		// quot = (newr / r)
		if ((r->leng == 1) && (r->nums[0] < 3))
		{
			bncopy(newr, quot);
			if (r->nums[0] > 1) { bnrshift(quot, 1); }
		}
		else { bndiv(newr, r, quot, temp); }
		// prev = s
		bncopy(s, prev);
		// s = (news - (quot * prev))
		bnzero(temp); bnmul(quot, prev, temp);
		bnsub(news, temp, s, 0);
		// news = prev
		bncopy(prev, news);
		// prev = r
		bncopy(r, prev);
		// r = (newr - (quot * prev))
		bnzero(temp); bnmul(quot, prev, temp);
		bnsub(newr, temp, r, 0);
		// newr = prev
		bncopy(prev, newr);
	}
	if (news->sign == 1)
	{
		// news = news + b
		bnadd(news, b, news, 0);
	}
	bncopy(news, g);
	bnfree(s); bnfree(news);
	bnfree(r); bnfree(newr);
	bnfree(prev); bnfree(quot); bnfree(temp);
}

// modular square root

int sqrtmod(bnum *a, bnum *p, bnum *r)
{
	bnum *o = bninit(1);
	o->nums[0] = 1; o->leng = 1; o->sign = 0;
	
	// legendre symbol
	// define if a is a quadratic residue modulo odd prime
	
	// g = (p - 1) / 2
	
	// p - 1
	bnum *qq = bndup(p);
	bnsub(qq, o, qq, 1);
	// (p - 1) / 2
	bnum *g = bndup(qq);
	bnrshift(g, 1);
	
	// l = pow(a, g, p)
	
	// pow(a, g, p)
	bnum *l = bninit(max(max(a->size, g->size), p->size) * 3);
	bnpowmod(a, g, p, l);
	if (bncmp(l, qq) == 0)
	{
		bnfree(o); bnfree(qq); bnfree(g); bnfree(l);
		return -1;
	}
	
	// factor p - 1 on the form q * (2 ^ s) (with Q odd)
	// q = p - 1; s = 0
	bnum *q = bndup(qq);
	bnum *s = bninit(p->size);
	while ((q->nums[0] % 2) == 0)
	{
		// s += 1; q /= 2
		bnadd(s, o, s, 1);
		bnrshift(q, 1);
	}
	
	// select a z which is a quadratic non resudue modulo p
	// z = 1
	bnum *z = bninit(p->size);
	z->nums[0] = 1; z->leng = 1; z->sign = 0;
	while (1)
	{
		// while (lsym(z, p) != -1)
		bnpowmod(z, g, p, l);
		if (bncmp(l, qq) == 0) { break; }
		// z += 1
		bnadd(z, o, z, 1);
	}
	// c = pow(z, q, p)
	bnum *c = bninit(max(max(z->size, q->size), p->size) * 3);
	bnpowmod(z, q, p, c);
	
	// search for a solution
	// f = ((q + 1) / 2)
	bnum *f = bndup(q);
	bnadd(f, o, f, 1); bnrshift(f, 1);
	// x = pow(a, f, p)
	bnpowmod(a, f, p, r);
	// t = pow(a, q, p)
	bnum *t = bninit(max(max(a->size, q->size), p->size) * 3);
	bnpowmod(a, q, p, t);
	// m = s
	bnum *m = bninit(p->size), *i = bninit(p->size), *e = bninit(p->size);
	bncopy(s, m);
	// u = 2
	bnum *u = bninit(1);
	u->nums[0] = 2; u->leng = 1; u->sign = 0;
	bnum *b = bninit(p->size * 4), *v = bninit(p->size * 4), *w = bninit(p->size * 4);
	while ((t->leng > 1) || (t->nums[0] != 1))
	{
		// find the lowest i such that t ^ (2 ^ i) = 1
		// i = 1; e = 2
		i->nums[0] = 1; i->leng = 1; i->sign = 0;
		e->nums[0] = 2; e->leng = 1; e->sign = 0;
		while (bncmp(i, m) < 0)
		{
			bnpowmod(t, e, p, l);
			if ((l->leng == 1) && (l->nums[0] == 1)) { break; }
			bnlshift(e, 1);
			bnadd(i, o, i, 1);
		}
		// update next value to iterate
		// (m - i - 1)
		bnsub(m, i, v, 0);
		bnsub(v, o, v, 0);
		// 2 ^ (m - i - 1)
		bnpowmod(u, v, p, l);
		// b = (c ^ (2 ^ (m - i - 1))) % p
		bnpowmod(c, l, p, b);
		// x = ((x * b) % p)
		bnzero(v); bnmul(r, b, v);
		bndiv(v, p, w, r);
		// b = (b * b) % p
		bnzero(v); bnmul(b, b, v);
		bndiv(v, p, w, b);
		// t = ((t * b) % p)
		bnzero(v); bnmul(t, b, v);
		bndiv(v, p, w, t);
		// c = b; m = i
		bncopy(b, c);
		bncopy(i, m);
	}
	
	bnfree(o); bnfree(qq); bnfree(g); bnfree(l);
	bnfree(q); bnfree(s); bnfree(z); bnfree(c);
	bnfree(f); bnfree(t); bnfree(m);
	bnfree(i); bnfree(e); bnfree(b);
	bnfree(u); bnfree(v); bnfree(w);
	
	// r = [x, p - x]
	
	return 0;
}

// montgomery curve arithmetic

void nmod(bnum *a, bnum *b)
{
	if (a->sign == 1) { bnadd(b, a, a, 0); }
}

void pdub(ecc *p, ecc *r, ect *t)
{
	// printf("2P=\n");
	
	// l = 3*x^2 + 2*a*x + 1 / 2*b*y
	
	// x^2
	bnzero(t->w); bnmul(p->x, p->x, t->w); bndiv(t->w, p->p, t->t, t->v);
	// 3*x^2
	bnadd(t->v, t->v, t->w, 1); bnadd(t->w, t->v, t->w, 1);
	// 2*a*x
	bnzero(t->h); bnmul(p->a, p->x, t->h); bnlshift(t->h, 1);
	bndiv(t->h, p->p, t->t, t->u); nmod(t->u, p->p);
	// 3*x^2 + 2*a*x + 1
	bnadd(t->w, t->u, t->g, 1);
	int x, o = 1;
	for (x = 0; (x < (t->g)->leng) && (o == 1); ++x)
	{
		o = 0; if ((t->g)->nums[x] == 0xffffffff) { o = 1; } (t->g)->nums[x] += 1;
	}
	if (o == 1) { (t->g)->nums[x] = 1; (t->g)->leng += 1; }
	bndiv(t->g, p->p, t->t, t->yr);
	// 1 / 2*b*y
	bnzero(t->w); bnmul(p->b, p->y, t->w); bnlshift(t->w, 1);
	bndiv(t->w, p->p, t->t, t->u); nmod(t->u, p->p); egcd(t->u, p->p, t->i);
	// 3*x^2 + 2*a*x + 1 / 2*b*y
	bnzero(t->w); bnmul(t->yr, t->i, t->w); bndiv(t->w, p->p, t->t, t->s);
	
	// xr = b*l^2 - a - 2*x
	
	// l^2
	bnzero(t->g); bnmul(t->s, t->s, t->g);
	// b*l^2 - a
	bnzero(t->w); bnmul(p->b, t->g, t->w);
	bnsub(t->w, p->a, t->w, 0);
	// 2*x
	bnzero(t->h); bncopy(p->x, t->h); bnlshift(t->h, 1);
	// b*l^2 - a - 2*x
	bnsub(t->w, t->h, t->w, 0);
	bndiv(t->w, p->p, t->t, t->xr); nmod(t->xr, p->p);
	
	// yr = ((3*x + a) * l) - b*l^3 - y
	
	// (3*x + a) * l
	bnadd(t->h, p->x, t->w, 1); bnadd(t->w, p->a, t->w, 0);
	bndiv(t->w, p->p, t->t, t->u); nmod(t->u, p->p);
	bnzero(t->w); bnmul(t->u, t->s, t->w);
	// l^3
	bncopy(t->g, t->h);
	bnzero(t->g); bnmul(t->h, t->s, t->g); bndiv(t->g, p->p, t->t, t->u);
	// b*l^3
	bnzero(t->h); bnmul(p->b, t->u, t->h);
	bndiv(t->h, p->p, t->t, t->u); nmod(t->u, p->p);
	// ((3*x + a) * l) - b*l^3 - y
	bnsub(t->w, t->u, t->w, 0); bnsub(t->w, p->y, t->w, 0);
	bndiv(t->w, p->p, t->t, t->yr); nmod(t->yr, p->p);
	
	(t->xr)->leng = (p->p)->size; bncopy(t->xr, r->x);
	(t->yr)->leng = (p->p)->size; bncopy(t->yr, r->y);
}

void padd(ecc *p, ecc *q, ecc *r, ect *t)
{
	// printf("P+Q=\n");
	
	// l = (Qy - Py) / (Qx - Px)
	
	// Qy - Py
	bnsub(q->y, p->y, t->yr, 0);
	// Qx - Px
	bnsub(q->x, p->x, t->xr, 0);
	bndiv(t->xr, p->p, t->t, t->u); nmod(t->u, p->p);
	// 1 / (Qx - Px)
	egcd(t->u, p->p, t->i);
	// (Qy - Py) / (Qx - Px)
	bnzero(t->w); bnmul(t->yr, t->i, t->w);
	bndiv(t->w, p->p, t->t, t->s);
	
	// xr = b*l^2 - a - Px - Qx
	
	// b*l^2 - a - Px - Qx
	bnzero(t->w); bnmul(t->s, t->s, t->w);
	bnzero(t->g); bnmul(p->b, t->w, t->g);
	bnsub(t->g, p->a, t->g, 0);
	bnsub(t->g, p->x, t->g, 0);
	bnsub(t->g, q->x, t->g, 0);
	bndiv(t->g, p->p, t->t, t->xr); nmod(t->xr, p->p);
	
	// yr = ((2*Px + Qx + a) * l) - b*l^3 - Py
	
	// 2*Px + Qx + a
	bnadd(p->x, p->x, t->t, 1);
	bnadd(t->t, q->x, t->t, 1);
	bnadd(t->t, p->a, t->t, 0);
	bndiv(t->t, p->p, t->u, t->v); nmod(t->v, p->p);
	// (2*Px + Qx + a) * l
	bnzero(t->w); bnmul(t->v, t->s, t->w); bndiv(t->w, p->p, t->t, t->u);
	// b*l^3
	bnzero(t->w); bnmul(t->s, t->s, t->w);
	bnzero(t->g); bnmul(t->w, t->s, t->g); bndiv(t->g, p->p, t->t, t->v);
	bnzero(t->w); bnmul(t->v, p->b, t->w);
	// ((2*Px + Qx + a) * l) - b*l^3 - Py
	bnsub(t->u, t->w, t->v, 0);
	bnsub(t->v, p->y, t->v, 0);
	bndiv(t->v, p->p, t->t, t->yr); nmod(t->yr, p->p);
	
	(t->xr)->leng = (p->p)->size; bncopy(t->xr, r->x);
	(t->yr)->leng = (p->p)->size; bncopy(t->yr, r->y);
}

void pmul(bnum *m, ecc *p, ecc *r)
{
	int init = 0;
	bnum *mul = bndup(m);
	ecc *b = ecdup(p);
	ect *t = etinit(p);
	while ((mul->leng > 1) || (mul->nums[0] > 0))
	{
		if ((mul->nums[0] % 2) == 1)
		{
			if (init == 0)
			{
				bnfree(r->x); r->x = bndup(b->x);
				bnfree(r->y); r->y = bndup(b->y);
			}
			else
			{
				padd(r, b, r, t);
			}
			init = 1;
		}
		pdub(b, b, t);
		bnrshift(mul, 1);
	}
	bnfree(mul);
	ecfree(b);
	etfree(t);
}