/*

*******************************************************************************

\file ecp.c

\brief Elliptic curves over prime fields

\project bee2 [cryptographic library]

\author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]

\created 2012.06.26

\version 2016.07.05

\license This program is released under the GNU General Public License 

version 3. See Copyright Notices in bee2/info.h.

*******************************************************************************

*/

#include "bee2/core/mem.h"

#include "bee2/core/stack.h"

#include "bee2/core/util.h"

#include "bee2/math/ecp.h"

#include "bee2/math/gfp.h"

#include "bee2/math/pri.h"

#include "bee2/math/ww.h"

#include "bee2/math/zz.h"

/*

*******************************************************************************

Общие положения

Ссылки на все реализованные алгоритмы имеются на сайте

http://www.hyperelliptic.org/efd. Там же можно найти соглашения по

обозначению сложности алгоритмов. В этих обозначениях фигурируют

следующие формальные выражения:

	add -- сложение или вычитание \mod p,

	c -- умножение на малую константу c (2,3,...) \mod p,

	half -- умножение на 2^{-1} \mod p,

	*A -- умножение на коэффициент A \mod p,

	S -- возведение в квадрат \mod p,

	M -- умножение \mod p,

	D -- деление \mod p.

При общей оценке сложности считается, что 1D = 100M и 1S = 1M.

Аддитивные операции игнорируются.

Используются обозначения:

	A <- A + A -- сложение аффинных точек,

	A <- 2A -- удвоение аффинной точки;

	P <- P + P -- сложение проективных точек;

	P <- P + A -- добавление к проективной точке аффинной;

	P <- 2P -- удвоение проективной точки;

	P <- 2A -- удвоение аффинной точки с переходом к проективным координатам.

*******************************************************************************

*/

#define ecpSeemsOnA(a, ec)\

	(zmIsIn(ecX(a), (ec)->f) && zmIsIn(ecY(a, (ec)->f->n), (ec)->f))

#define ecpSeemsOn3(a, ec)\

	(ecpSeemsOnA(a, ec) && zmIsIn(ecZ(a, (ec)->f->n), (ec)->f))

/*

*******************************************************************************

Якобиановы координаты:

	x = X / Z^2, y = Y / Z^3,

	-(X : Y : Z) = (X : -Y : Z).

В функции ecpDblJ() выполняется удвоение P <- 2P. Реализован

алгоритм dbl-1998-hnm [Hasegawa T., Nakajima J., Matsui M. A Practical

Implementation of Elliptic Curve Cryptosystems over GF(p) on a

16-bit Microcomputer. Public Key Cryptography, 1998].

Сложность алгоритма:

	3M + 6S + 1*A + 1half + 6add + 3*2 \approx 9M.

\todo Сравнить dbl-1998-hnm с dbl-2007-bl, сложность которого

	1M + 8S + 1A + 10add + 1*8 + 2*2 + 1*3.

В функции ecpDblJA3() выполняется удвоение P <- 2P для случая A = -3.

Реализован алгоритм dbl-1998-hnm2 [Hasegawa T., Nakajima J., Matsui M.

A Practical Implementation of Elliptic Curve Cryptosystems over GF(p)

on a 16-bit Microcomputer. Public Key Cryptography, 1998].

Сложность алгоритма:

	4M + 4S + 1*half + 7add + 3*2 \approx 8M.

\todo Сравнить dbl-1998-hnm2 с dbl-2001-b, сложность которого

	3M +5S + 8add + 1*4 + 2*8 + 1*3.

В функции ecpDblAJ() выполняется удвоение P <- 2A. Реализован алгоритм 

mdbl-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма:

	1M + 5S + 7add + 1*8 + 3*2 + 1*3 \approx 6M.

В функции ecpAddJ() выполняется сложение P <- P + P. Реализован алгоритм 

add-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма:

	11M + 5S + 9add + 4*2 \approx 16M.

\todo Сравнить add-2007-bl с add-1998-cmo-2, сложность которого

	12M + 4S + 6add + 1*2.

В функции ecpAddAJ() выполняется сложение P <- P + A.

Реализован алгоритм madd-2004-hmv [Hankerson D., Menezes A., Vanstone S.

Guide to Elliptic Curve Cryptography, Springer, 2004].

Сложность алгоритма:

	8M + 3S + 6add + 1*2 \approx 11M.

В функции ecpTplJ() выполняется утроение P <- 3P. Реализован алгоритм 

tpl-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма

	5M + 10S + 1*A + 15add + 2*4 + 1*6 + 1*8 + 1*16 + 1*3 \approx 15M.

В функции ecpTplJA3() выполняется утроение P <- 3P для случая A = -3. 

Реализован алгоритм tpl-2007-bl-2 [Bernstein-Lange, 2007]. 

Сложность алгоритма

	7M + 7S + 13add + 2*4 + 1*8 + 1*12 + 1*16 + 1*3 \approx 14M.

Целевые функции ci(l), определенные в описании реализации ecMul() в ec.c, 

принимают следующий вид:

	c1(l) = l/3 11;

	c2(l, w) = 103 + (2^{w-2} - 2)102 + l/(w + 1) 11;

	c3(l, w) = 6 + (2^{w-2} - 2)16 + l/(w + 1) 16.

Расчеты показывают, что

	с1(l) <= min_w c3(l), l <= 81,

	min_w c3(l, w) <= m	in_w c2(l, w), l <= 899.

Поэтому для практически используемых размерностей l (192 <= l <= 571)

первые две стратегии являются проигрышными. Реализована только стратегия 3.

\todo Сравнить madd-2004-hmv с madd-2007-bl, сложность которого

	7M + 4S + 9add + 1*4 + 3*2.

\todo При выполнении функции ecpAddAJ() может выясниться, что

складываются одинаковые точки. В этом случае вызывается функция

ecpDblAJ(). Подумать, как "помочь" этой функции, передав уже

рассчитанные значения, например, za^2. Аналогичное замечание касается

функции ecpAddJ().

*******************************************************************************

*/

// [3n]b <- [2n]a (P <- A)

{

	// pre

	// xb <- xa

	// yb <- ya

	// zb <- unity

}

// [2n]b <- [3n]a (A <- P)

{

	// переменные в stack

	// pre

	// a == O => b <- O

	// t1 <- za^{-1}

	// t2 <- t1^2

	// xb <- xa t2

	// t2 <- t1 t2

	// yb <- ya t2

	// b != O

}

{

}

// [3n]b <- -[3n]a (P <- -P)

{

	// pre

	// xb <- xa

	// yb <- -ya

	// zb <- za

}

// [3n]b <- 2[3n]a (P <- 2P)

{

	// переменные в stack

	// pre

	// za == 0 или ya == 0? => b <- O

	{

	}

	// t1 <- za^2

	// zb <- ya za

	// zb <- 2 zb

	// t1 <- t1^2

	// t1 <- A t1

	// t2 <- xa^2

	// t1 <- t1 + t2

	// t2 <- 2 t2

	// t1 <- t1 + t2

	// yb <- 2 ya

	// yb <- yb^2

	// t2 <- yb^2

	// t2 <- t2 / 2

	// yb <- yb xa

	// xb <- t1^2

	// xb <- xb - yb

	// xb <- xb - yb

	// yb <- yb - xb

	// yb <- yb t1

	// yb <- yb - t2

}

{

}

// [3n]b <- 2[3n]a (P <- 2P, A = -3)

{

	// переменные в stack

	// pre

	// za == 0 или ya == 0? => b <- O

	{

	}

	// t1 <- za^2

	// zb <- ya za

	// zb <- 2 zb

	// t2 <- xa - t1

	// t1 <- xa + t1

	// t2 <- t1 t2

	// t1 <- 2 t2

	// t1 <- t1 + t2

	// yb <- 2 ya

	// yb <- yb^2

	// t2 <- yb^2

	// t2 <- t2 / 2

	// yb <- yb xa

	// xb <- t1^2

	// xb <- xb - yb

	// xb <- xb - yb

	// yb <- yb - xb

	// yb <- yb t1

	// yb <- yb - t2

}

{

}

// [3n]b <- 2[2n]a (P <- 2A)

{

	// переменные в stack

	// pre

	// ya == 0? => b <- O

	{

	}

	// t1 <- xa^2 [X1^2 = XX]

	// t2 <- ya^2 [Y1^2 = YY]

	// t3 <- t2^2 [YY^2 = YYYY]

	// t2 <- t2 + xa [X1 + YY]

	// t2 <- t2^2 [(X1 + YY)^2]

	// t2 <- t2 - t1 [(X1 + YY)^2 - XX]

	// t2 <- t2 - t3 [(X1 + YY)^2 - XX - YYYY]

	// t2 <- 2 t2 [2((X1 + YY)^2 - XX - YYYY) = S]

	// t4 <- 2 t1 [2 XX]

	// t4 <- t4 + t1 [3 XX]

	// t4 <- t4 + A [3 XX + A = M]

	// t1 <- 2 t2 [2S]

	// xb <- t4^2 [M^2]

	// xb <- xb - t1 [M^2 - 2S = T]

	// zb <- 2 ya [2Y1]

	// t2 <- t2 - xb [S - T]

	// yb <- t4 t2 [M(S - T)]

	// t3 <- 2 t3 [2 YYYY]

	// t3 <- 2 t3 [4 YYYY]

	// t3 <- 2 t3 [8 YYYY]

	// yb <- yb - t3 [M(S - T) - 8 YYYY]

}

{

}

// [3n]c <- [3n]a + [3n]b (P <- P + P)

	void* stack)

{

	// переменные в stack

	// pre

	// a == O => c <- b

	{

	}

	// b == O => c <- a

	{

	}

	// t1 <- za^2 [Z1Z1]

	// t2 <- zb^2 [Z2Z2]

	// t3 <- zb t2 [Z2 Z2Z2]

	// t3 <- ya t3 [Y1 Z2 Z2Z2 = S1]

	// t4 <- za t1 [Z1 Z1Z1]

	// t4 <- yb t4 [Y2 Z1 Z1Z1 = S2]

	// zc <- za + zb [Z1 + Z2]

	// zc <- zc^2 [(Z1 + Z2)^2]

	// zc <- zc - t1 [(Z1 + Z2)^2 - Z1Z1]

	// zc <- zc - t2 [(Z1 + Z2)^2 - Z1Z1 - Z2Z2]

	// t1 <- xb t1 [X1 Z2Z2 = U2]

	// t2 <- xa t2 [X2 Z1Z1 = U1]

	// t1 <- t1 - t2 [U2 - U1 = H]

	// t1 == 0 => xa zb^2 == xb za^2

	{

		// t3 == t4 => ya zb^3 == yb za^3 => a == b => c <- 2a

		// t3 != t4 => a == -b => c <- O

		else

	}

	// zc <- zc t1 [((Z1 + Z2)^2 - Z1Z1 - Z2Z2)H = Z3]

	// t4 <- t4 - t3 [S2 - S1]

	// t4 <- 2 t4 [2(S2 - S1) = r]

	// yc <- 2 t1 [2H]

	// yc <- yc^2 [(2H)^2 = I]

	// t1 <- t1 yc [H I = J]

	// yc <- t2 yc [U1 I = V]

	// t2 <- 2 yc [2 V]

	// xc <- t4^2 [r^2]

	// xc <- xc - t1 [r^2 - J]

	// xc <- xc - t2 [r^2 - J - 2V = X3]

	// yc <- yc - xc [V - X3]

	// yc <- t4 yc [r(V - X3)]

	// t3 <- 2 t3 [2S1]

	// t3 <- t3 t1 [2S1 J]

	// yc <- yc - t3 [r(V - X3) - 2 S1 J]

}

{

			f_deep,

			ecpDblJ_deep(n, f_deep));

}

// [3n]c <- [3n]a + [2n]b (P <- P + A)

	void* stack)

{

	// переменные в stack

	// pre

	// a == O => c <- (xb : yb : 1)

	{

	}

	// t1 <- za^2

	// t2 <- t1 za

	// t1 <- t1 xb

	// t2 <- t2 yb

	// t1 <- t1 - xa

	// t2 <- t2 - ya

	// t1 == 0?

	{

		// t2 == 0 => c <- 2(xb : yb : 1)

		// t2 != 0 => c <- O

		else

	}

	// zc <- t1 za

	// t3 <- t1^2

	// t4 <- t1 t3

	// t3 <- t3 xa

	// t1 <- 2 t3

	// xc <- t2^2

	// xc <- xc - t1

	// xc <- xc - t4

	// t3 <- t3 - xc

	// t3 <- t3 t2

	// t4 <- t4 ya

	// yc <- t3 - t4

}

{

			f_deep,

			ecpDblAJ_deep(n, f_deep));

}

// [3n]c <- [3n]a - [3n]b (P <- P - P)

	void* stack)

{

	// переменные в stack

	// pre

	// t <- -b

	// c <- a + t

}

{

}

// [3n]c <- [3n]a - [2n]b (P <- P - A)

	void* stack)

{

	// переменные в stack

	// pre

	// t <- -b

	// c <- a + t