/*

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

\file ec2.c

\brief Elliptic curves over binary fields

\project bee2 [cryptographic library]

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

\created 2012.06.26

\version 2014.07.15

\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/ec2.h"

#include "bee2/math/gf2.h"

#include "bee2/math/pri.h"

#include "bee2/math/ww.h"

#include "bee2/math/zz.h"

/*

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

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

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

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

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

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

	add -- сложение или вычитание в GF(2^m),

	c -- умножение на малую константу c в GF(2^m),

	*A -- умножение на коэффициент A в GF(2^m),

	*B -- умножение на коэффициент B в GF(2^m),

	S -- возведение в квадрат в GF(2^m),

	M -- умножение в GF(2^m),

	D -- деление в GF(2^m).

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

Аддитивные операции игнорируются. В общем случае 1*A = 1M, но в наиболее

распространенных на практике случаях A \in {0, 1} и 1*A = 0M.

\warning Соотношение 1D = 24M получено экспериментальным путем

на платформе x86 и возможно требует пересмотра.

На упомянутом сайте http://www.hyperelliptic.org/efd используется другое

соотношение: 1D = 10M. Считаем его слишком оптимистичным

(по отношению к D).

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

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

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

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

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

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

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

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

*/

#define ec2SeemsOnA(a, ec)\

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

#define ec2SeemsOn3(a, ec)\

	(ec2SeemsOnA(a, ec) && gf2IsIn(ecZ(a, (ec)->f->n), (ec)->f))

/*

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

Кривая в проективных координатах Лопеса -- Дахаба (LD):

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

	O = (1 : 0 : 0),

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

\warning Ошибка в книге [Hankerson D., Menezes A., Vanstone S. Guide to

Elliptic Curve Cryptography, Springer, 2004] при определении обратной

точки в LD-координатах.

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

dbl-2005-l [Lange, 2005]. Сложность алгоритма:

	4M + 4S + 1*A + 5add \approx 5M,

причем умножение на A не выполняется, если A \in {0, 1}.

\todo Ссравнить с алгоритмом dbl-2005-dl.

В функции ec2DblALD() выполняется удвоение P <- 2A (Z-координата

проективной точки равняется 1). Реализован алгоритм

mdbl-2005-dl [Doche-Lange, 2005]. Сложность алгоритма:

	1M + 3S + 1*A + 1*B + 4add \approx 3M.

причем умножение на A не выполняется, если A \in {0, 1}.

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

Реализован алгоритм add-2005-dl [Doche-Lange-Takagi, 2005].

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

	13M + 4S + 9add \approx 13M.

В функции ec2AddALD() выполняется сложение P <- P + A (Z-координата

второго слагаемого равняется 1).

Реализован алгоритм madd-2005-dl [Doche, Lange, Al-Daoude, 2005].

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

		8M + 5S + 1*A + 9add \approx 9M,

причем умножение на A не выполняется, если A \in {0, 1}.

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

принимают следующий вид (считаем, что коэффициент A \in {0, 1}):

	c1(l) = l/3 8;

	c2(l, w) = 26 + (2^{w-2} - 2)26 + l/(w + 1) 8;

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

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

	с1(l) <= min_w c2(l), l <= 39,

	min_w c2(l, w) <= min_w c3(l, w).

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

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

\todo Реализовать быстрые формулы для особенных B:

B = 1 (кривые Коблица),	известен \sqrt{B}.

\todo Реализовать редакции функций с A \in {0, 1}.

\todo Исследовать расширенные LD-координаты (дополнительно поддерживается Z^2).

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

*/

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

	void* stack)

{

	// pre

	// xb <- xa

	// yb <- ya

	// zb <- 1

}

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

{

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

	// pre

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

	// t1 <- za^{-1}

	// xb <- xa t1

	// t1 <- t1^2

	// yb <- ya t1

	// b != O

}

{

}

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

{

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

	// pre

	// t1 <- xa * za

	// b <- (xa, ya + t1, za)

}

{

}

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

{

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

	// pre

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

	{

	}

	// t1 <- xa za [A]

	// zb <- t1^2 [A^2]

	// t2 <- xa^2 [B]

	// xb <- ya + t2 [C]

	// t1 <- t1 xb [D]

	// xb <- xb^2 + t1 [C^2 + D]

	// t1 <- t1 + zb [Z3 + D]

	// yb <- t2^2 zb [B^2 Z3]

	// xb <- xb + A * zb [C^2 + D + a2 * Z3]

	{

	}

	// t1 <- t1 xb [(Z3 + D)X3]

	// yb <- yb + t1 [(Z3 + D)X3 + B^2 Z3]

}

{

}

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

{

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

	// pre

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

	{

	}

	// zb <- xa^2 [C]

	// xb <- zb^2 + B [C^2 + a6]

	// yb <- ya^2 + B [Y1^2 + a6]

	// yb <- yb + A zb [Y1^2 + a2*Z3 + a6]

	{

	}

	// yb <- yb xb [(Y1^2 + a2*Z3 + a6) * X3]

	// t1 <- B zb [a6 * Z3]

	// yb <- yb + t1 [(Y1^2 + a2*Z3 + a6) * X3 + a6 * Z3]

}

{

}

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

	const ec_o* ec, void* stack)

{

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

	// pre

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

	{

	}

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

	{

	}

	// t1 <- xa zb [A]

	// t2 <- xb za [B]

	// t3 <- ya zb^2 [G]

	// t4 <- yb za^2 [H]

	// A == B => a == \pm b

	{

		// t3 == t4 => a == b => c <- 2a

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

		else

	}

	// t5 <- t1 + t2 [E]

	// t6 <- t3 + t4 [I]

	// t5 <- t5 t6 [J]

	// xc <- t1^2 [C]

	// yc <- t2^2 [D]

	// t6 <- xc + yc [ec->f]

	// zc <- t6 za zb [ec->f * Z1 * Z2]

	// t4 <- t1 (t4 + yc) [A * (H + D)]

	// xc <- t2 (xc + t3) + t4 [B * (C + G) + A * (H + D)]

	// t1 <- t1 t5 [A * J]

	// t3 <- t3 t6 [ec->f * G]

	// t1 <- (t1 + t3) t6 [(A * J + ec->f * G) * ec->f]

	// yc <- (t5 + zc) xc [(J + Z3) * X3]

	// yc <- yc + t1 [(A * J + ec->f * G) * ec->f + (J + Z3) * X3]

}

{

			f_deep,

			ec2DblLD_deep(n, f_deep));

}

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

	const ec_o* ec, void* stack)

{

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

	// pre

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

	{

	}

	// t1 <- ya + yb za^2 [A]

	// t2 <- xa + xb za [B]

	// t2 == 0 => a == \pm b

	{

		// t1 == 0 => a == b => c <- 2b

		// t1 != 0 => a == -b => c <- O

		else

	}

	// t3 <- t2 za [C]

	// zc <- t3^2 [C^2]

	// t4 <- xb zc [D]

	// yc <- xb + yb [X2 + Y2]

	// xc <- t2^2 + t1 + A t3 [B^2 + A + a2 * C]

	{

	}

	// xc <- xc t3 + t1^2 [C * (A + B^2 + a2 * C) + A^2]

	// yc <- yc zc^2 [(Y2 + X2) * Z3^2]

	// t4 <- t4 + xc [D + X3]

	// t1 <- t1 t3 + zc [A * C + Z3]

	// t1 <- t1 t4 [(D + X3)(A * C + Z3)]

	// yc <- yc + t1 [(D + X3)(A * C + Z3) + (Y2 + X2) * Z3^2]

}

{

			f_deep,

			ec2DblALD_deep(n, f_deep));

}

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

	const ec_o* ec, void* stack)

{

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

	// pre

	// t <- -b

	// c <- a + t

}

{

			f_deep,

			ec2AddLD_deep(n, f_deep));

}

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

	const ec_o* ec, void* stack)

{

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

	// pre

	// t <- -b

	// c <- a + t

}

{

}

	void* stack)

{

	// обнулить

	// зафикисровать размерности

	// запомнить базовое поле

	// сохранить коэффициенты

	// подготовить буферы для описания группы точек

	// настроить интерфейсы

		ec2ToALD_deep(f->n, f->deep),

		ec2NegLD_deep(f->n, f->deep),

		ec2AddLD_deep(f->n, f->deep),

		ec2AddALD_deep(f->n, f->deep),

		ec2SubLD_deep(f->n, f->deep),

		ec2SubALD_deep(f->n, f->deep),

		ec2DblLD_deep(f->n, f->deep),

		ec2DblALD_deep(f->n, f->deep));

	// настроить заголовок

	// все нормально

}

{

}

{

		ec2ToALD_deep(n, f_deep),

		ec2NegLD_deep(n, f_deep),

		ec2AddLD_deep(n, f_deep),

		ec2AddALD_deep(n, f_deep),

		ec2SubLD_deep(n, f_deep),

		ec2SubALD_deep(n, f_deep),

		ec2DblLD_deep(n, f_deep),

		ec2DblALD_deep(n, f_deep));

}

/*

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

Свойства кривой

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

*/

{

	// кривая работоспособна? поле ec->f корректно?

	// ec->deep >= ec->f->deep?

	// A, B \in ec->f, B != 0?

}

{

}

{

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