椭圆曲线加密

csuzhucong

已于 2023-10-18 23:01:39 修改

阅读量5.4k

点赞数 7

文章标签：人工智能算法

于 2021-04-12 20:55:12 首次发布

本文链接：https://blog.csdn.net/nameofcsdn/article/details/115627882

版权

3，Weierstrass方程（魏尔斯特拉斯方程）

五，Montgomery曲线、X25519

1，Montgomery曲线

（1）方程

（2）Weierstrass和Montgomery互转

（2）Montgomery和Edwards互转

〇，知识背景

本文涉及到的知识背景：

（1）三次曲线、三次方程的判别式

（2）群、环、域、阿贝尔群、循环群、有限域、生成子群、同构群、群阶、拉格朗日定理

（3）同余、完全剩余系、逆元、加法生成元、欧拉φ函数、平方剩余、欧拉准则、二次互反律

（4）埃及乘法、离散对数问题、公钥体系

一，椭圆

椭圆方程： $\left ( \frac{x}{a} \right )^2+\left ( \frac{y}{b} \right )^2=1$ ，其中a>b>0

椭圆面积：S=PI * a * b

椭圆周长比较复杂，没有初等公式。

椭圆的周长可以用完全椭圆积分表示，其中就会涉及到一个形如 $\sqrt{t^3+...}$ 的式子。

二，椭圆曲线

1，一般椭圆曲线

（1）方程

方程 $y^2+uxy+vy=ax^3+bx^2+cx+d$ 所描述的曲线，如果处处光滑（没有尖点、孤点、自交），那么就是椭圆曲线

（2）非椭圆曲线

2个反例：尖点和自交

下面的Montgomery曲线一章中，有孤点的例子。

（3）规约成常用椭圆曲线

$y^2+uxy+vy=ax^3+bx^2+cx+d$ 可以变成 $(y+\frac{u}{2}x+\frac{v}{2})=ax^3+(b+\frac{u^2}{4})x^2+(c+\frac{uv}{2})x+d+\frac{v^2}{4}$

（4）退化为圆锥曲线

2，常用椭圆曲线

$y^2=ax^3+bx^2+cx+d$ 即 $y=\pm \sqrt{ax^3+bx^2+cx+d}$ ，其中a不为0

非退化的椭圆曲线都可以规约成这个形式。

因为 $\sqrt{ax^3+bx^2+cx+d}$ 的形式和完全椭圆积分中的式子一致，所以被称为椭圆曲线，实际上和椭圆的关系并不大。

后面只讨论常用椭圆曲线，加密算法中应该也只涉及这个。

右边为三次项：

3，Weierstrass方程（魏尔斯特拉斯方程）

常用椭圆曲线都可以规约成如下的简单形式：

$y^2=x^3+ax+b$

这个好像是叫Weierstrass normal form吧，没太仔细查，不确定。

实际上，三次方程最常用的求根公式，也是先规约成这个形式，然后根据韦达定理可以推出求根公式。

4，椭圆曲线的判别式

以Weierstrass方程为例，下文也默认都讨论Weierstrass方程

判别式为 $\Delta = 4a^3+27b^2$ ，和三次方程的判别式相同，参考求根公式

这其实很好理解，椭圆曲线的形状，就是取决于它的对称轴和它的交点情况。

如果判别式为0，那么要么有三重根（尖点），要么有两重根（自交）

如果判别式不为0，那么就一定是光滑的，也就是椭圆曲线了。

5，椭圆曲线上的加法

以 $y^2=x^3+ax+b$ 为例

我们来构造一个加法运算，让加法符合群的规则数学与泛型编程（2）群论_加群所有有限阶元素构成g的一个子群-CSDN博客

过曲线上的两点A、B画一条直线，找到直线与椭圆曲线的交点，交点关于x轴对称位置的点，定义为A+B，即为加法。

感谢图源椭圆曲线加密算法 - 简书

作直线关于x轴（即椭圆曲线的对称轴）的对称直线

于是A+B+B'=A

为了让加法符合群，我们使用结合律，得到B+B'为单位元素0

于是，我们把B的对称点定义为-B

如果A、B重合，那么就过点A作切线。

6，加法群分析

OK，有了这个基本运算之后，我们再来全面分析，这样的加法能不能构成一个群？

（1）直线和曲线的交点

如果直线是垂直于x轴的，那么直线和曲线没有其他交点，也没有相切属性。

如果直线不是垂直于x轴，那么有方程组

$\left\{\begin{matrix} y^2=ax^3+bx^2+cx+d\\ y=kx+t\: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \end{matrix}\right.$

既然已经有2个交点了，就一定会有第三个交点，但是可能会有重合的解，即切线的情况。

（2）封闭性

要想任意两个点相加之和仍然还在集合内，我们除了椭圆曲线上的所有点之外，还增加一个无穷远点

即，任意垂直于x轴的直线，与曲线的交点都在无穷远点，那么显然，它就是单位元素。

对于在x轴上的那些点D（至少1个，可能有多个），过点D作切线，也相交于无穷远点，

也就是说，D=-D，注意，这并不能推导出D=0，

类似的群，由平面的旋转方式所形成的群：群元素是指平面围绕原点的各种旋转方式，群运算是对两种旋转方式所做的合成，单位元素是旋转0度，D是旋转180度，D+D就是旋转0度

（3）群的运算和常量

群运算就是如上的加法运算（A+A就是作切线），逆运算就是取对称点（x轴上的点和无穷远点的对称点都是自身）,

唯一的单位元素就是无穷远点

（4）群的公理

单位元素公理和消去公理已经很明显了，主要是结合性公理。

A+B+C=A+(B+C)

结合性的证明很困难，中文搜了几十个网页都找不到任何关于如何证明的内容，连思路都没有任何博主提到，甚至抄来抄去之后很多博主以为这个是显然的，其实这个是很不显然的。

7，椭圆曲线加法的结合性证明

我们使用Weierstrass normal form，它和直线相交的方程组是

$\left\{\begin{matrix} y^2=x^3+ax+b\\ y=kx+b\: \: \: \: \: \: \: \! \: \end{matrix}\right.$

根据韦达定理，三个根的和为 $-k^{2}$

由此，如果已知直线和曲线相交的两个点，就可以轻松求出第三个点的坐标。

虽然群不需要交换性，只有阿贝尔群需要交换性，但是椭圆曲线加法的交换性太显然了，所以可以直接用。

下面证明结合性associativity：

其中是曲线上的加法运算。

设P(x1,y1) Q(x2,y2) R(x3,y3)

设P+Q的坐标是（s,t），那么

其中Det是行列式determinant

设P+Q+R的坐标是（CC，DD），那么

PS：要证明结合性，只需要证明，CC和DD对于x1 y1、x2 y2、x3 y3是轮换对称的。

设Q+R的坐标是（u，v），那么

设Q+R+P的坐标是（FF，GG），那么

至此，我们已经得到了（CC，DD）和（FF，GG），他们是轮换式，即把x1 y1和x3 y3互换，

要证明结合性，只需要证明CC-FF=0，DD-GG=0

只需要证明最后一行这个式子为0即可，然而暴力展开的话大概有十几万项，吓人猪心！

除以

利用matlab计算除的结果发现可以整除，而这个式子是恒为0的，PQ不是恒为0的，所以就证明出CC-FF=0，同理可得DD-GG=0

代码：

syms a b c d e f
A=f*(b-a)+e*(c-a)+d*(b-c);
B=f*(b-a)-e*(c-a)+d*(b-c);
P=(d-e)^2-(a-b)^2*(a+b+c);
Q=(f-e)^2-(b-c)^2*(a+b+c);
S=((a-b)*A*Q)^2-((c-b)*B*P)^2+2*P*Q*((e-d)*A*Q-(e-f)*B*P)+(P*Q)^2*2*(a-c);
T=(b-c)*d^2+(c-a)*e^2+(a-b)*f^2+(a-b)*(b-c)*(c-a)*(a+b+c);
[r,q] = polynomialReduce(S,T)

运行：

r = 0
q =

- 2*a^6*b^3 + 6*a^6*b^2*c - 6*a^6*b*c^2 + 2*a^6*c^3 - 6*a^5*b^3*c + 18*a^5*b^2*c^2 - 18*a^5*b*c^3 + 3*a^5*b*e^2 - 6*a^5*b*e*f + 3*a^5*b*f^2 + 6*a^5*c^4 - 3*a^5*c*e^2 + 6*a^5*c*e*f - 3*a^5*c*f^2 + 6*a^4*b^5 - 12*a^4*b^4*c - 6*a^4*b^3*c^2 + 30*a^4*b^2*c^3 - 2*a^4*b^2*e^2 + 8*a^4*b^2*e*f - 6*a^4*b^2*f^2 - 24*a^4*b*c^4 + 4*a^4*b*c*e^2 - 16*a^4*b*c*e*f + 12*a^4*b*c*f^2 + 6*a^4*c^5 - 2*a^4*c^2*e^2 + 8*a^4*c^2*e*f - 6*a^4*c^2*f^2 + 12*a^3*b^5*c - 24*a^3*b^4*c^2 - 2*a^3*b^3*c^3 + 3*a^3*b^3*d^2 - 8*a^3*b^3*d*e + 6*a^3*b^3*d*f - 3*a^3*b^3*e^2 + 2*a^3*b^3*e*f + 30*a^3*b^2*c^4 - 9*a^3*b^2*c*d^2 + 18*a^3*b^2*c*d*e - 12*a^3*b^2*c*d*f - 4*a^3*b^2*c*e^2 + 22*a^3*b^2*c*e*f - 15*a^3*b^2*c*f^2 - 18*a^3*b*c^5 + 9*a^3*b*c^2*d^2 - 12*a^3*b*c^2*d*e + 6*a^3*b*c^2*d*f + 5*a^3*b*c^2*e^2 - 26*a^3*b*c^2*e*f + 18*a^3*b*c^2*f^2 + 2*a^3*c^6 - 3*a^3*c^3*d^2 + 2*a^3*c^3*d*e + 2*a^3*c^3*e^2 + 2*a^3*c^3*e*f - 3*a^3*c^3*f^2 + 2*a^3*d*e^3 - 6*a^3*d*e^2*f + 6*a^3*d*e*f^2 - 2*a^3*d*f^3 - 3*a^3*e^4 + 8*a^3*e^3*f - 6*a^3*e^2*f^2 + a^3*f^4 - 6*a^2*b^7 + 6*a^2*b^6*c + 18*a^2*b^5*c^2 - 24*a^2*b^4*c^3 + 2*a^2*b^4*d*e - 6*a^2*b^4*d*f + 4*a^2*b^4*e^2 - 6*a^2*b^4*e*f + 6*a^2*b^4*f^2 - 6*a^2*b^3*c^4 + 6*a^2*b^3*c*d^2 - 14*a^2*b^3*c*d*e + 18*a^2*b^3*c*d*f - a^2*b^3*c*e^2 - 12*a^2*b^3*c*e*f + 3*a^2*b^3*c*f^2 + 18*a^2*b^2*c^5 - 18*a^2*b^2*c^2*d^2 + 30*a^2*b^2*c^2*d*e - 18*a^2*b^2*c^2*d*f - 6*a^2*b^2*c^2*e^2 + 30*a^2*b^2*c^2*e*f - 18*a^2*b^2*c^2*f^2 - 6*a^2*b*c^6 + 18*a^2*b*c^3*d^2 - 26*a^2*b*c^3*d*e + 6*a^2*b*c^3*d*f + 5*a^2*b*c^3*e^2 - 12*a^2*b*c^3*e*f + 9*a^2*b*c^3*f^2 - 3*a^2*b*d^2*e^2 + 6*a^2*b*d^2*e*f - 3*a^2*b*d^2*f^2 + 4*a^2*b*d*e^3 - 8*a^2*b*d*e^2*f + 4*a^2*b*d*e*f^2 + 2*a^2*b*e^4 - 4*a^2*b*e^3*f - a^2*b*e^2*f^2 + 6*a^2*b*e*f^3 - 3*a^2*b*f^4 - 6*a^2*c^4*d^2 + 8*a^2*c^4*d*e - 2*a^2*c^4*e^2 + 3*a^2*c*d^2*e^2 - 6*a^2*c*d^2*e*f + 3*a^2*c*d^2*f^2 - 4*a^2*c*d*e^3 + 8*a^2*c*d*e^2*f - 4*a^2*c*d*e*f^2 + a^2*c*e^4 - 2*a^2*c*e^3*f + a^2*c*e^2*f^2 - 6*a*b^7*c + 6*a*b^6*c^2 + 12*a*b^5*c^3 - 3*a*b^5*d^2 + 8*a*b^5*d*e - 6*a*b^5*d*f + 4*a*b^5*e*f - 3*a*b^5*f^2 - 12*a*b^4*c^4 + 6*a*b^4*c*d^2 - 8*a*b^4*c*d*e + 4*a*b^4*c*e^2 - 8*a*b^4*c*e*f + 6*a*b^4*c*f^2 - 6*a*b^3*c^5 + 3*a*b^3*c^2*d^2 - 12*a*b^3*c^2*d*e + 18*a*b^3*c^2*d*f - a*b^3*c^2*e^2 - 14*a*b^3*c^2*e*f + 6*a*b^3*c^2*f^2 + 6*a*b^2*c^6 - 15*a*b^2*c^3*d^2 + 22*a*b^2*c^3*d*e - 12*a*b^2*c^3*d*f - 4*a*b^2*c^3*e^2 + 18*a*b^2*c^3*e*f - 9*a*b^2*c^3*f^2 + a*b^2*d^2*e^2 - 4*a*b^2*d^2*e*f + 3*a*b^2*d^2*f^2 - 4*a*b^2*d*e^3 + 18*a*b^2*d*e^2*f - 20*a*b^2*d*e*f^2 + 6*a*b^2*d*f^3 - 8*a*b^2*e^3*f + 17*a*b^2*e^2*f^2 - 12*a*b^2*e*f^3 + 3*a*b^2*f^4 + 12*a*b*c^4*d^2 - 16*a*b*c^4*d*e + 4*a*b*c^4*e^2 - 2*a*b*c*d^2*e^2 + 8*a*b*c*d^2*e*f - 6*a*b*c*d^2*f^2 - 8*a*b*c*d*e^2*f + 8*a*b*c*d*e*f^2 + 2*a*b*c*e^4 - 2*a*b*c*e^2*f^2 - 3*a*c^5*d^2 + 6*a*c^5*d*e - 3*a*c^5*e^2 + a*c^2*d^2*e^2 - 4*a*c^2*d^2*e*f + 3*a*c^2*d^2*f^2 - 2*a*c^2*d*e^3 + 8*a*c^2*d*e^2*f - 6*a*c^2*d*e*f^2 + a*c^2*e^4 - 4*a*c^2*e^3*f + 3*a*c^2*e^2*f^2 + 2*b^9 - 6*b^7*c^2 - 2*b^6*d*e + 6*b^6*d*f - 2*b^6*e^2 - 2*b^6*e*f + 6*b^5*c^4 - 3*b^5*c*d^2 + 4*b^5*c*d*e - 6*b^5*c*d*f + 8*b^5*c*e*f - 3*b^5*c*f^2 + 6*b^4*c^2*d^2 - 6*b^4*c^2*d*e - 6*b^4*c^2*d*f + 4*b^4*c^2*e^2 + 2*b^4*c^2*e*f - 2*b^3*c^6 + 2*b^3*c^3*d*e + 6*b^3*c^3*d*f - 3*b^3*c^3*e^2 - 8*b^3*c^3*e*f + 3*b^3*c^3*f^2 - b^3*d^4 + 6*b^3*d^3*e - 4*b^3*d^3*f - 10*b^3*d^2*e^2 + 10*b^3*d^2*e*f + 8*b^3*d*e^3 - 16*b^3*d*e^2*f + 10*b^3*d*e*f^2 - 4*b^3*d*f^3 - 2*b^3*e^4 + 8*b^3*e^3*f - 10*b^3*e^2*f^2 + 6*b^3*e*f^3 - b^3*f^4 - 6*b^2*c^4*d^2 + 8*b^2*c^4*d*e - 2*b^2*c^4*e^2 + 3*b^2*c*d^4 - 12*b^2*c*d^3*e + 6*b^2*c*d^3*f + 17*b^2*c*d^2*e^2 - 20*b^2*c*d^2*e*f + 3*b^2*c*d^2*f^2 - 8*b^2*c*d*e^3 + 18*b^2*c*d*e^2*f - 4*b^2*c*d*e*f^2 - 4*b^2*c*e^3*f + b^2*c*e^2*f^2 + 3*b*c^5*d^2 - 6*b*c^5*d*e + 3*b*c^5*e^2 - 3*b*c^2*d^4 + 6*b*c^2*d^3*e - b*c^2*d^2*e^2 + 4*b*c^2*d^2*e*f - 3*b*c^2*d^2*f^2 - 4*b*c^2*d*e^3 - 8*b*c^2*d*e^2*f + 6*b*c^2*d*e*f^2 + 2*b*c^2*e^4 + 4*b*c^2*e^3*f - 3*b*c^2*e^2*f^2 + c^3*d^4 - 2*c^3*d^3*f - 6*c^3*d^2*e^2 + 6*c^3*d^2*e*f + 8*c^3*d*e^3 - 6*c^3*d*e^2*f - 3*c^3*e^4 + 2*c^3*e^3*f - 2*d^3*e^3 + 6*d^3*e^2*f - 6*d^3*e*f^2 + 2*d^3*f^3 + 6*d^2*e^4 - 18*d^2*e^3*f + 18*d^2*e^2*f^2 - 6*d^2*e*f^3 - 6*d*e^5 + 18*d*e^4*f - 18*d*e^3*f^2 + 6*d*e^2*f^3 + 2*e^6 - 6*e^5*f + 6*e^4*f^2 - 2*e^3*f^3

8，阿贝尔群

证明出结合性之后，椭圆曲线上的加法就是群了，由于与生俱来的交换性，使得它是一个阿贝尔群。

三，素数域上的椭圆曲线

1，有限域、素数域

域

素数域Fp：由p个元素（从0到p-1）和数论加减乘除构成的域

2，素数域上的椭圆曲线

椭圆曲线Ep(a,b)， $y^2\equiv x^3+ax+b(modp)$ ， x,y∈[0,p-1]

其中p为质数，非负整数a、b小于p且满足 $4a^3+27b^2\neq 0(modp)$

素数域上的椭圆曲线，除了满足同余方程的所有点之外，同样需要额外加一个无穷远点，即单位元素O点

写个简单的程序，来感受一下椭圆曲线上的点的分布。


#include<iostream>
#include<windows.h>
#include<map>
using namespace std;

#define CIN(x) while (!(cin >> x)) { \
        cin.clear();      \
        cin.ignore();     \
    }
#define CIN2(x, y) CIN(x)CIN(y)
#define CIN3(x, y, z) CIN(x)CIN(y)CIN(z)

const int M = 10000;

bool check(int p)
{
    if(p<3 || p>M)return false;
    if(p%2==0)return false;
    for(int i=3;i*i<=p;i+=2)if(p%i==0)return false;
    return true;
}

void play()
{
    cout<<"y^2=x^3+ax+b(mod p)，输入p,a,b\n";
    int p,a,b;
    CIN3(p,a,b);
    if(!check(p))return;
    a=(a%p+p)%p,b=(b%p+p)%p;
    map<int,int>modSqrt;
    for(int i=0;i<=p/2;i++)modSqrt[i*i%p]=i+1;
    for(int i=0;i<p;i++){
        int x=(i*i*i+a*i+b)%p;
        if(modSqrt[x]==0)continue;
        if(modSqrt[x]==1)cout<<i<<" "<<0<<"    ";
        else cout<<i<<" "<<modSqrt[x]-1<<"    "<<i<<" "<<p-modSqrt[x]+1<<"    " ;
    }
}

int main() {
    for(int i=0;i<10;i++){
        play();
        sleep(5);
    }
    return 0;
}

输入23 1 3，输出：

0 7 0 16 2 6 2 17 4 5 4 18 5 8 5 15 6 8 6 15 7 10 7 13 10 1 10 22 12 8 12 15 14 1 14 22 15 9 15 14 19 2 19 21 21 4 21 19 22 1 22 22

再加上O点，一共27个点。

改成可视化版本：

void play()
{
    cout<<"y^2=x^3+ax+b(mod p)，输入p,a,b\n";
    int p,a,b;
    CIN3(p,a,b);
    if(!check(p))return;
    a=(a%p+p)%p,b=(b%p+p)%p;
    map<int,int>modSqrt;
    for(int i=0;i<=p/2;i++)modSqrt[i*i%p]=i+1;
    if(p>30)return;
    int pla[30][30]={0};
    for(int i=0;i<p;i++){
        int x=(i*i*i+a*i+b)%p;
        if(modSqrt[x]==0)continue;
        pla[i][modSqrt[x]-1]=1,pla[i][p-modSqrt[x]+1]=1;
    }
    for(int j=p-1;j>=0;j--){
        for(int i=0;i<p;i++){
            if(pla[i][j])cout<<'*';
            else cout<<'.';
        }
        cout<<endl;
    }
}

把输出结果的截图修饰一下得到：

上面13个点和下面13个点是对称的，O点单独存在。

准确的说，对于这个p行p列的正方形来说，上面的p-1行是上下对称的。

换句话说，如果把上面的（p-1）/ 2 行循环挪到下面，就会得到以y=0这一行为对称轴的，一个上下对称的p行p列的正方形。

3，素数域上椭圆曲线的加法

（1）定义交点

和实数域上一样，P和Q连成的直线和曲线相交于R，P(x1,y1), Q(x2,y2), R(x3,y3)三个点都在直线y=kx+b上，

那么x1+x2+x3=k^2，其中k是斜率。

（2）定义逆元素

同样的，互为上下对称的2个点A和A'互为相反元素，A+A'=O

（3）定义加法

若P(x1,y1), Q(x2,y2), R(x3,y3)满足P+Q=R，

则 $x_3=k^2-x_1-x_2,\, y_3=k(x_1-x_3)-y_1$ ，其中k是直线PQ的斜率，

（3.1）若 $P\neq Q$ ，则 $k=\frac{y_2-y_1}{x_2-x_1}$

若x1=x2，则k就是无穷，R是无穷远点O

（3.2）若P=Q，则PQ是切线， $k=\frac{3x_1^2+a}{2y_1}$

若y1=0，则k就是无穷，R是无穷远点O

（4）阿贝尔群

素数域上的加法，构成的也是阿贝尔群，而且是循环群。

（5）demo程序

#include<iostream>
#include<map>
using namespace std;

long long get_mi(long long n,int k,int p)
{
	if (k == 0)return 1;
	long long r = get_mi(n, k / 2,p) % p;
	r = (r*r) % p;
	if (k % 2)r = (r*n) % p;
	return r;
}

long long reverse(int a,int p)
{
    return get_mi(a,p-2,p);
}

#define CIN(x) while (!(cin >> x)) { \
        cin.clear();      \
        cin.ignore();     \
    }
#define CIN2(x, y) CIN(x)CIN(y)
#define CIN3(x, y, z) CIN(x)CIN(y)CIN(z)

const int M = 10000;

bool check(int p)
{
    if(p<3 || p>M)return false;
    if(p%2==0)return false;
    for(int i=3;i*i<=p;i+=2)if(p%i==0)return false;
    return true;
}

bool isO(int x,int y)
{
    return x==0 && y==0;
}

void setO(int &x,int &y)
{
    x=0,y=0;
}

int p,a,b;

void play()
{
    cout<<"y^2=x^3+ax+b(mod p),input p,a,b\n";

    CIN3(p,a,b);
    if(!check(p))return;
    a=(a%p+p)%p,b=(b%p+p)%p;
    map<int,int>modSqrt;
    for(int i=0;i<=p/2;i++)modSqrt[i*i%p]=i+1;
    for(int i=0;i<p;i++){
        int x=(i*i*i+a*i+b)%p;
        if(modSqrt[x]==0)continue;
        if(modSqrt[x]==1)cout<<i<<" "<<0<<"    ";
        else cout<<i<<" "<<modSqrt[x]-1<<"    "<<i<<" "<<p-modSqrt[x]+1<<"    " ;
    }
}

void add(int x1,int y1,int x2,int y2,int &x3,int &y3)
{
    int k;
    if(isO(x1,y1))x3=x2,y3=y2;
    else if(isO(x2,y2))x3=x1,y3=y1;
    else if(x1==x2&&y1==y2){
        if(y1==0)setO(x3,y3);
        else k=(x1*x1*3+a)%p*(p+1)/2*reverse(y1,p)%p,x3=((k*k-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }else{
        if(x1==x2)setO(x3,y3);
        else k=(y2-y1)*reverse(x2-x1,p)%p,x3=((k*k-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }
}

void play2(int n)
{
    int x1,y1,x2,y2,x3,y3,k;
    while(n--){
        cout<<"\n input x1,y1,x2,y2\n";
        CIN2(x1,y1);
        CIN2(x2,y2);
        add(x1,y1,x2,y2,x3,y3);
        cout<<x3<<" "<<y3;
    }
}

int main()
{
    play();
    play2(10);
	return 0;
}

运行：

其中，无穷远点可以有很多表示方法，为了方便可以在有限域内选择不在曲线上的任意一点，比如当b不为0时，（0,0）就可以用来表示无穷远点。

4，素数域上点的个数——Hasse定理

corollary 2 是 theorem 1在k=1时的特例的一个推论

5，生成子群

加法构成的整个集合是有限群，而由单个元素的数乘构成的是它的子群。

根据拉格朗日定理数学与泛型编程（2）群论_加群所有有限阶元素构成g的一个子群-CSDN博客

子群的阶一定是群阶的约数。

还是以E(23,1,3)为例，我们看下每个元素的生成子群的阶：

#include<iostream>
#include<map>
#include<vector>
using namespace std;

long long get_mi(long long n,int k,int p)
{
	if (k == 0)return 1;
	long long r = get_mi(n, k / 2,p) % p;
	r = (r*r) % p;
	if (k % 2)r = (r*n) % p;
	return r;
}

long long reverse(int a,int p)
{
    return get_mi(a,p-2,p);
}

#define CIN(x) while (!(cin >> x)) { \
        cin.clear();      \
        cin.ignore();     \
    }
#define CIN2(x, y) CIN(x)CIN(y)
#define CIN3(x, y, z) CIN(x)CIN(y)CIN(z)

const int M = 10000;

bool check(int p)
{
    if(p<3 || p>M)return false;
    if(p%2==0)return false;
    for(int i=3;i*i<=p;i+=2)if(p%i==0)return false;
    return true;
}

bool isO(int x,int y)
{
    return x==0 && y==0;
}

void setO(int &x,int &y)
{
    x=0,y=0;
}

int p,a,b;
vector<int>vx;
vector<int>vy;


void play()
{
    cout<<"y^2=x^3+ax+b(mod p),input p,a,b\n";
    CIN3(p,a,b);
    if(!check(p))return;
    a=(a%p+p)%p,b=(b%p+p)%p;
    map<int,int>modSqrt;
    for(int i=0;i<=p/2;i++)modSqrt[i*i%p]=i+1;
    for(int i=0;i<p;i++){
        int x=(i*i*i+a*i+b)%p;
        if(modSqrt[x]==0)continue;
        if(modSqrt[x]==1)vx.push_back(i),vy.push_back(0);
        else vx.push_back(i),vx.push_back(i),vy.push_back(modSqrt[x]-1),vy.push_back(p-modSqrt[x]+1);
    }
    int x,y;
    setO(x,y);
    vx.push_back(x),vy.push_back(y);
}

void add(int x1,int y1,int x2,int y2,int &x3,int &y3)
{
    int k;
    if(isO(x1,y1))x3=x2,y3=y2;
    else if(isO(x2,y2))x3=x1,y3=y1;
    else if(x1==x2&&y1==y2){
        if(y1==0)setO(x3,y3);
        else k=(x1*x1*3+a)%p*(p+1)/2*reverse(y1,p)%p,x3=((k*k-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }else{
        if(x1==x2)setO(x3,y3);
        else k=(y2-y1)*reverse(x2-x1,p)%p,x3=((k*k-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }
}

void getJ()
{
    int x1,y1,x2,y2;
    for(int i=0;i<vx.size();i++) {
        x1=x2=vx[i],y1=y2=vy[i];
        int s=1;
        while(!isO(x1,y1)){
            s++;
            add(x1,y1,x2,y2,x1,y1);
        }
        cout<<x2<<" "<<y2<<"   "<<s<<"      ";
    }
}

int main()
{
    play();
    getJ();
	return 0;
}

输入23 1 3 ，输出：

0 7 27 0 16 27 2 6 27 2 17 27 4 5 9 4 18 9 5 8 27 5 15 27 6 8 27
6 15 27 7 10 27 7 13 27 10 1 9 10 22 9 12 8 3 12 15 3 14 1 27
14 22 27 15 9 27 15 14 27 19 2 9 19 21 9 21 4 27 21 19 27 22 1 27
22 22 27 0 0 1

不难发现规律，所有元素的阶都是大群阶的约数，我们把大群阶的约数除以元素的阶叫做该元素的协因子co-factor

在加密的时候，一般选择协因子最小的元素。

对于E(23,1,3)，对于每个约数，阶为这个约数的元素的数目是欧拉函数，即有φ(27)个元素的阶是27

对于其他曲线，不一定就是循环群。

PS： $\sum _{d|n}\varphi (d)=n$

四，ECC椭圆曲线加密

1，数学计算

因为有限域上的加法是阿贝尔群，所以数乘的计算可以采用埃及乘法数学与泛型编程（1）埃及乘法、加法链_数学与泛型编程 pdf-CSDN博客

所以，对于任意点P，任意正整数k，计算Q=kP是很容易的，而根据P和Q计算k是很难的，目前的最好算法是O（sqrt(p)）的时间复杂度，BSGS算法离散对数求解、BSGS算法_离散对数问题_csuzhucong的博客-CSDN博客

2，加密方案

P、Q作为公钥公开，k作为私钥保留

发送方把信息转化为曲线上的一个点A，并产生一个随机整数x，并计算出2个点：xP和xQ

再把B=xP和C=A+xQ这2个点发送给接收方。

接收方收到B和C之后，根据C-kB=A+xQ-kxP=A即可算出A

3，破解方案

（1）根据P和Q算出私钥k，对使用该公钥进行加密的信息一次性全部破解

（2）根据B=xP算出随机加密种子x，只能破解用x进行加密的信息，即只能破解单独一份信息。

现在一般p取二进制200位，sqrt(p)就是100位，即10^30

现在的超级计算机的算力大概10^17，破解需要30万年

4，constant-time设计

接收方完成计算C-kB的时间只和k有关，

准确的说，快速幂算法的时间只和k的位数有关，所以破解方可以根据算法运行时间缩小范围，从而破解k，即时间旁路攻击

为了避免时间旁路攻击，可以手动控制循环次数，即把k理解为固定长度的数，如X25519中无论k的实际长度多长，都视为256位。

五，Montgomery曲线、X25519

1，Montgomery曲线

（1）方程

曲线方程：By^2=x(x^2+Ax+1)，其中B(A^2-4)不为0

Montgomery曲线有三种形态：

当A^2-4=0时，曲线会有孤点（1,0）或者（-1,0）

Montgomery曲线应该不存在尖点或者自交的形态。

（2）Weierstrass和Montgomery互转

虽然x(x^2+Ax+1)=0这个三次方程可以直接求解，不需要求根公式，

但是我们还是可以用求根公式的规约法，令x = u - A/3，把它化成关于u的没有二次项的三次多项式，

然后令By^2=v^2，这样就把Montgomery化成了Weierstrass，反之同理。

（3）加法

若P(x1,y1), Q(x2,y2), R(x3,y3)满足P+Q=R，

则 $x_3=Bk^2-A-x_1-x_2,\, y_3=k(x_1-x_3)-y_1$ ，其中k是直线PQ的斜率，

（3.1）若 $P\neq Q$ ，则 $k=\frac{y_2-y_1}{x_2-x_1}$

若x1=x2，则k就是无穷，R是无穷远点O

（3.2）若P=Q，则PQ是切线， $k=\frac{3x_1^2+2Ax_1+1}{2By_1}$

若y1=0，则k就是无穷，R是无穷远点O

2，有限域上的Montgomery曲线

（1）唯一固定点

（0,0）唯一能确定在曲线上的点，而这个点的阶显然是2

（2）阶

In any finite field, B(A+ 2), B(A−2), and A^2 −4 cannot all be nonsquares simultaneously

论文：https://eprint.iacr.org/2017/212.pdf

因为B(A+ 2)*B(A−2)=B^2*（A^2 −4）所以三者至少有1个是平方剩余

（2.1）如果B(A+2)是平方剩余，那么有限域上必有一点X的横坐标是1，而过X的切线恰好经过（0,0），所以X的阶就是4

（2.2）如果B(A-2)是平方剩余，那么有限域上必有一点X的横坐标是-1，而过X的切线恰好经过（0,0），所以X的阶就是4

（2.3）如果B(A+2)和B(A-2)都不是平方剩余，则A^2 −4是平方剩余，方程x^2+Ax+1就有2个不同的解，而且都不为0，即有3个不同的点都满足纵坐标为0

那么这3个点加上无穷远点可以构成一个子群，这个子群的阶是4

所以大群阶一定是4的倍数，即曲线上所有点的个数是4的倍数0

我的代码：

#include<iostream>
#include<map>
#include<vector>
using namespace std;

// 省略代码同上

#define XOO -12345678
bool isO(int x,int y)
{
    return x==XOO && y==XOO;
}

void setO(int &x,int &y)
{
    x=XOO,y=XOO;
}

int p,a,b;
vector<int>vx;
vector<int>vy;


void play()
{
    cout<<"By^2=x^3+Ax^2+x(mod p),input p,A,B\n";
    CIN3(p,a,b);
    if(!check(p))return;
    a=(a%p+p)%p,b=(b%p+p)%p;
    map<int,int>modSqrt;
    for(int i=0;i<=p/2;i++)modSqrt[i*i%p]=i+1;
    for(int i=0;i<p;i++){
        int x=(i*i*i+a*i*i+i)%p;
        x=x*reverse(b,p)%p;
        if(modSqrt[x]==0)continue;
        if(modSqrt[x]==1)vx.push_back(i),vy.push_back(0);
        else vx.push_back(i),vx.push_back(i),vy.push_back(modSqrt[x]-1),vy.push_back(p-modSqrt[x]+1);
    }
    int x,y;
    setO(x,y);
    vx.push_back(x),vy.push_back(y);
}

void add(int x1,int y1,int x2,int y2,int &x3,int &y3)
{
    int k;
    if(isO(x1,y1))x3=x2,y3=y2;
    else if(isO(x2,y2))x3=x1,y3=y1;
    else if(x1==x2&&y1==y2){
        if(y1==0)setO(x3,y3);
        else k=(x1*x1*3+x1*a*2+1)%p*(p+1)/2*reverse(y1,p)%p*reverse(b,p)%p,x3=((k*k*b-a-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }else{
        if(x1==x2)setO(x3,y3);
        else k=(y2-y1)*reverse(x2-x1,p)%p,x3=((k*k*b-a-x1-x2)%p+p)%p,y3=((k*(x1-x3)-y1)%p+p)%p;
    }
}

void getJ()
{
    int x1,y1,x2,y2;
    cout<<vx.size()<<endl;
    for(int i=0;i<vx.size();i++) {
        x1=x2=vx[i],y1=y2=vy[i];
        int s=1;
        while(!isO(x1,y1)){
            s++;
            add(x1,y1,x2,y2,x1,y1);
        }
        cout<<x2<<" "<<y2<<"   "<<s<<"      ";
    }
}

void play2(int n)
{
    int x1,y1,x2,y2,x3,y3;
    while(n--){
        cout<<"\n input x1,y1,x2,y2\n";
        CIN2(x1,y1);
        CIN2(x2,y2);
        add(x1,y1,x2,y2,x3,y3);
        cout<<x3<<" "<<y3;
    }
}

int main()
{
    play();
    getJ();
    play2(40);
	return 0;
}

代码基本上和Weierstrass 方程的代码一样，只需要把参数改一改即可，另外无穷远点也要换个表示法。

运行：

By^2=x^3+Ax^2+x(mod p),input p,A,B
23 11 5
32
0 0   2      3 9   16      3 14   16      4 9   16      4 14   16      5 9   8      5 14   8      6 2   16      6 21   16
7 7   4      7 16   4      8 1   16      8 22   16      9 6   16      9 17   16      10 10   4      10 13   4
11 1   8      11 22   8      13 8   16      13 15   16      14 7   8      14 16   8      15 0   2      16 1   16
16 22   16      18 11   16      18 12   16      20 0   2      21 4   8      21 19   8      -12345678 -12345678   1

3，X25519（Curve25519）

（1）方程

X25519的方程：y^2=x^3+486662x^2+x，p=2^255-19

rfc7748：RFC 7748 - Elliptic Curves for Security

下面是宏观图和微观图

486662这个值应该是基于2^255-19算出来的，寻找的目标就是不能太大（算法要快）而且阶要有大因子

因为阶J一定是4的倍数，所以阶的因子最多可能是J/4，其次最多可能是J/8，而J的大小只能在极其有限的范围内波动（Hasse定理）

所以最理想的方案是J=4p，其中p是素数，p约等于2^253，或者J=8p，其中p是素数，p约等于2^252

486662对应的是后者这个方案。

（2）基点

基点取x=9，y=14781619447589544791020593568409986887264606134616475288964881837755586237401

阶是2^252+27742317777372353535851937790883648493（素数），协因子co-factor=8

可以用python命令行验证一下基点坐标值：

六，Edwards曲线、ED25519

1，Edwards曲线

（1）方程

$au^2 + v^2 = 1 + du^2v^2$

可以规约成 $u^2 + v^2 = 1 + du^2v^2$ （如果a是平方剩余）

或 $-u^2 + v^2 = 1 + du^2v^2$ （如果-a是平方剩余）

（2）Montgomery和Edwards互转

可以规约成 $u^2 + v^2 = 1 + du^2v^2,d=(A-2)/(A+2)$ （如果B(A+2)是平方剩余）

或 $-u^2 + v^2 = 1 + du^2v^2,d=-(A-2)/(A+2)$ （如果-B(A+2)是平方剩余）

2，ED25519

（1）方程

直接Montgomery中的把X25519转化成Edwards：

y^2=x^3+486662x^2+x，p=2^255-19 即A=486662，B=1

所以 a=486664, d=486660

即 $486664u^2 + v^2 = 1 + 486660u^2v^2$

能不能规约取决于486664是不是平方剩余

$\left ( \frac{486664}{p} \right )=\left ( \frac{2}{p} \right )\left ( \frac{127}{p} \right )\left ( \frac{479}{p} \right )$

根据二次互反律 $\left(\frac{2}{{p}}\right)=(-1)^{\frac{p^{2}-1}{8}}$ ， $\left(\frac{{p}}{q}\right)=(-1)^{\frac{(p-1)(q-1)}{4}}\left(\frac{q}{p}\right)$

可以得到 $\left ( \frac{2}{p} \right )=-1, \left ( \frac{127}{p} \right )=\left ( \frac{p}{127} \right )=\left ( \frac{116}{127} \right ), \left ( \frac{479}{p} \right )=\left ( \frac{p}{479} \right )=\left ( \frac{373}{479} \right )$

所以 $\left ( \frac{486664}{p} \right )=-\left ( \frac{116}{127} \right )\left ( \frac{373}{479} \right )=-\left ( \frac{127}{29} \right )\left ( \frac{479}{373} \right )=-1*(-1)*1=1$