ECC非对称加密算法

椭圆曲线

椭圆曲线在代数上的表示是下面这个方程：
y² = x³ + ax + b
其中，a = 0, b = 7 (比特币系统所使用的版本)，它的图形如下：

椭圆曲线有一些很有用的特征

• 一条非垂直的直线与椭圆曲线相交于两点，若这两点均不是切点，则曲线上必有第三点与那条直线相交
• 过曲线上任意一点的非垂直切线与该曲线必有且仅有另一个交点。

利用这些特征，我们可以定义两种运算：“异点相加”和“同点加倍”。

“异点相加”, P + Q = r, 定义为：r为r’基于x轴的反射点（对称点）。其中，R’为包含P和Q的直线与曲线的第三个交点,如图上所示。

同样，“同点加倍”，P + P = r, 定义为：作一条过P点的切线，先求出该切线与曲线的另一交点R’,再计算r‘基于x轴的反射点r。

求r 坐标,得到一个非常美的结果

当p!=q

当p=q

有限域

同时，并不是所有的椭圆曲线都适合加密。y²=x³+ax+b是一类可以用来加密的椭圆曲线，也是最为简单的一类。下面我们就把y²=x³+ax+b 这条曲线定义在Fp（模p剩余类构成的域）上：
选择两个满足下列条件的小于p(p为素数)的非负整数a、b
4a³+27b²≠0　(mod p)
则满足下列方程的所有点(x,y)，再加上 无穷远点O∞ ，构成一条椭圆曲线。
y²=x³+ax+b (mod p)
其中 x,y属于0到p-1间的整数，准确的说是模素数p剩余类，并将这条椭圆曲线记为Ep(a,b)。

这里为什么要加上无穷远点呢，无穷远点来自于射影平面，射影平面比欧式平面多了无穷远点，所有无穷远点构成无穷远直线，射影平面有一个重要假设：

平行线在无穷远处相较于一个点，即无穷远点O点，并且每条直线只有一个无穷远点

在椭圆曲线E_p(a,b)中p₁+r₁=O,p₁+O=p₁,p₂+r₂=O,p₂+O=p₂
所有椭圆线点按照P+Q=r算法构成加群
O为单元零元，p₁,r₁互为逆元，p₂,r₂互为逆元。

举个例子

令p = 71,a=0,b=7，曲线点已经离散了，但还是对称的，对称点互为逆元
加群有72个元素（加一个无穷远点）每个元素阶如下。

令p = 79,a=0,b=7，加群元素个数67（素数），素数阶群，每个元素的阶（除了单位元）都是67，都是群的生成元，计算出来结果

算法原理

考虑如下等式：K=kG [其中 K,G为Ep(a,b)上的点，k为小于n（n是点G的阶）的整数]，不难发现，给定k和G，根据加法法则，计算K很容易；但给定K和G，求k就相对困难了。这就是椭圆曲线加密算法采用的难题，我们把点G称为基点（base point）。

加解密流程：
1、用户A选定一条椭圆曲线Ep(a,b)，并取椭圆曲线上一点，作为基点G。
2、用户A选择一个私有密钥k，并生成公开密钥K=kG。
3、用户A将Ep(a,b)和点K，G传给用户B。
4、用户B接到信息后 ，将待传输的明文编码到Ep(a,b)上一点M（编码方法很多，这里不作讨论），并产生一个随机整数r。
5、用户B计算点C1=M+rK；C2=rG。
6、用户B将C1、C2传给用户A。
7、用户A接到信息后，计算C1-kC2，结果就是点M。因为 C1-kC2=M+rK-k(rG)=M+rK-r(kG)=M,
再对点M进行解码就可以得到明文。

G = (1,18)
prikey = 40
pubkey = get_add(G, prikey)
r = 16
M = (34,24)
C1 = get_r(M,get_add(pubkey, r))
C2 = get_add(G, r)
temp = get_add(C2,prikey)
get_r(C1,(temp[0], MOD-temp[1]))]))

最后一步解密结果和原文一样

[34, 24]

签名过程如下：

1、选择一条椭圆曲线Ep(a,b)，和基点G；
2、选择私有密钥k（k<n，n为G的阶），利用基点G计算公开密钥K=kG；
3、产生一个随机整数r（r<n），计算点R=rG；
4、将原数据和点R的坐标值x,y作为参数，计算SHA1做为hash，即Hash=SHA1(原数据,x,y)；
5、计算s≡r - Hash * k (mod n)
6、r和s做为签名值，如果r和s其中一个为0，重新从第3步开始执行

验证过程如下：
1、接受方在收到消息(m)和签名值(r,s)后，进行以下运算
2、计算：sG+H(m)K=(x1,y1), r1≡ x1 mod p。
3、验证等式：r1 ≡ r mod p。
4、如果等式成立，接受签名，否则签名无效。

import numpy as np
import  sys  
import matplotlib.pyplot as plt
def ecc_equation(a, b):
    def ecc(x):
        return x**3+a*x+7
    return ecc

# 求平方根
def get_square_root_mod(mod):
    def get_square_root(n):
        lr = []
        if n >= mod:
            n = n%mod
        for el in xrange(0, mod):
            if el**2%mod == n:
                return el
    return get_square_root
# y
def get_y(x):
    return [x,get_square_root(ecc(x))]

# 检查点是否在曲线上
def check_point(p, mod=MOD):
    if p[1]**2%mod == (ecc(p[0]))%mod:
        return True
    else:
        return False
# 求element在模mod剩余类群逆元
def invert(element, mod):
    if element >= mod:
        element = element%mod
    if element == 0:
        return None
    for index in xrange(1, mod):
        if element*index%mod == 1:
            return index
# 给出p,q求r=p+q
# if p != q
# c = (py-qy)/(px-qx)
# rx = c^2 - px-qx
# ry = c(px-rx)-py
# if p ==q
# c = (3px^2+a)/2py,rx = c^2-2px,ry=c(px-rx)-py
def get_r(p, q,mod=MOD, a=A):
    p = map(lambda x: x % mod, p)
    q = map(lambda x: x % mod, q)
    if p[0] == q[0] and (p[1]+q[1])%mod==0:
        #互为逆元点和为无穷远点，方便处理 记为[np.infty,np.infty]
        return [np.infty,np.infty]
    if p != q:
        c = (p[1]-q[1])*invert(p[0]-q[0], mod)%mod
    else:
        c = (3*p[0]**2+a)*invert(2*p[1],mod)%mod
    rx = (c**2-p[0]-q[0])%mod
    ry = (c*(p[0]-rx)-p[1])%mod
    return [rx,ry]

# 构建循环群
def cycle_group(generate_el):
    power = [generate_el]
    lr = generate_el
    while True:
        lr = get_r(generate_el, lr)
        power.append(lr)
        if lr == [np.infty, np.infty]:
            return power
# 求倍点
def get_add(G, multiple):
    lr = G
    for index in xrange(1, multiple):
        lr = get_r(lr, G)
    return lr
# 功能同上，利用同点加倍，性能更高
def get_multiple(G, multiple):
    if multiple%2 == 0:
        temp = get_multiple(G, multiple/2)
        return get_r(temp,temp)
    else:
        if multiple > 1:
            temp = get_multiple(G, multiple-1)
            return get_r(G, temp)
        else:
            return G
import numpy as np
import  sys  
import matplotlib.pyplot as plt

A = 0
B = 7
MOD = 79
#79
ecc = ecc_equation(A,B)
get_square_root = get_square_root_mod(MOD)
x = xrange(0,MOD)
y = map(get_y, x)
y = np.array(y)
total = y[np.where(y[:,1]>-1)]
total_inverse =map(lambda x:[x[0], MOD-x[1]],filter(filtery , total))
total = np.concatenate((total, total_inverse),axis=0)
#print total

plt.scatter(total[:,0], total[:,1])
plt.show()
print len(total)
print total[np.where(total[:,1]==0)]
print map(lambda x:len(x), map(cycle_group, total))