这篇博客分享了作者在毕设中实现的一些数论和椭圆曲线算法，包括有限域求逆元、快速幂、Legendre符号计算、Tonelli-Shanks开根、椭圆曲线类等。还介绍了numpy复数类的重载和有理点类，并提供了最大公因数、最小公倍数的计算方法。此外，博主探讨了穷举法生成椭圆曲线、类数算法以及Hilbert类域多项式生成等高级主题。
序言
毕设答辩算是告一段落了，大学生涯预告着结束了。
但这并不意味我将会闲暇下来，我还要继续考研，所以这篇文章编写可能不会花太多的时间。
我打算把我毕设中的一些算法发出来，供大家参考学习。
希望能帮到一些人。
当然中间的一些算法可能会有bug（我把比较稳定的算法都放在前面，一些不是很好用的我放到了后面来说。ps：因为不是所有算法都在毕设中用到了）
其实能够发现，前面的代码会比较规范，后面的代码就比较随性了。。。
算法目录
- 有限域求逆（function）
- 有限域快速幂（function）
- Legendre符号计算（function）
- Tonelli-Shanks有限域开根（function）
- 椭圆曲线类（class）
- numpy复数类的重载（class）
- 有理点类（class）
- 最大公因数
- 最小公倍数
- 数组最大公因数
- 数组最小公倍数
- 穷举法生成特定点数的椭圆曲线
- 类数算法
- Hilbert 类域多项式的生成
- j不变量的数值计算
- 曲线生成（CM法）
- 基于CRT算法的Hilbert类域多项式计算并生成曲线
- 有限域上多项式计算
有限域求逆元
def euc_inv(M: list) -> list:
'''a single step for inv_fin'''
x = M[0][0]
y = M[1][0]
if x == 0 or y == 0:
return M
if M[0][0] < 0 and M[1][0] < 0:
M[0] = [-M[0][0], -M[0][1], -M[0][2]]
M[1] = [-M[1][0], -M[1][1], -M[1][2]]
if abs(x) >= abs(y):
k = x // y
r = x - k * y
M[0] = [M[0][0] - k * M[1][0], M[0][1] - k * M[1][1], M[0][2] - k * M[1][2]]
return euc_inv(M)
else:
k = y // x
r = y - k * x
M[1] = [M[1][0] - k * M[0][0], M[1][1] - k * M[0][1], M[1][2] - k * M[0][2]]
return euc_inv(M)
def inv_fin(a: int, char: int) -> False or int:
'''get the inverse in finite field, Euclidean'''
if a % char == 0:
return False
if a % char == 1:
return 1
elif a % char == char - 1:
return -1
M = [[a, 1, 0], [char, 0, 1]]
M = euc_inv(M)
if abs(M[0][0]) == 1:
r = M[0][1]
return r
elif abs(M[1][0]) == 1:
r = M[1][1]
return r
else:
return False
有限域快速幂
def quick_powmod(a: int, b: int, c: int) -> int:
a = a % c
ans = 1
while b != 0:
if b & 1:
ans = (ans * a) % c
a = (a * a) % c
b >>= 1
return ans
Legendre符号计算
def legendre(n: int, char: int) -> 1 or -1 or None:
leg = quick_powmod(n, int((char - 1) // 2), char)
leg = leg % char
if n == 0 :
return 0
if leg == char - 1:
return -1
elif leg == 1:
return 1
else:
return NoneSolovay-Strassen大素数判定
def is_prime(n: int, r: bool = False, rng: int = 15) -> bool:
'''Solovay-Strassen Algorithm :to discriminate whether a number is prime'''
m = min(rng, int((n - 1) / 2))
a = [legendre(i, n) for i in range(1, m)]
if n <= 20:
if n in [2, 3, 5, 7, 11, 13, 17, 19]:
return True
else:
return False
else:
if None in a:
if r:
print(a)
return False
else:
if r:
print(a)
return TrueTonelli-Shanks有限域开根
def decom(p: int) -> (int, int):
'''for quadratic residue'''
Q, S = 1, 1
k = bin(p-1)
while True:
if k[-S] == str(1):
break
S += 1
S -= 1
Q = (p - 1) // 2 ** S
return Q, S
def quad_res(n: int, p: int) -> False or (int, int):
"""Tonelli-Shanks algorithm"""
if is_prime(p, rng=30) is False:
return False
chooseN = lambda p: [i for i in range(2, int(p / 2) + 1) if legendre(i, p) == -1][0]
if n == 0:
return 0, 0
if legendre(n, p) == -1:
return False
Q, S = decom(p)
if S == 1:
r = quick_powmod(n, int((p + 1) / 4), p)
return r, p - r
z = chooseN(p)
C = quick_powmod(z, Q, p)
R = quick_powmod(n, (Q + 1) // 2, p)
t = quick_powmod(n, Q, p)
M = S
while True:
if t % p == 1:
return R, p - R
else:
i = [item for item in range(M) if
quick_powmod(t, quick_powmod(2, item, p), p) % p == 1]
i = i[0]
b = quick_powmod(C, quick_powmod(2, M - i - 1, p), p)
R = R * b % p
t = t * b ** 2 % p
C = b ** 2 % p
M = i % p椭圆曲线类
该类中包含以下函数
- 给横坐标求纵坐标
- 给横坐标返回有理点
- 明文嵌入（3和4可能会有一些问题）
- 明文恢复
- 生成可除多项式
- 有理点阶数计算
- 有理点个数计算（穷举法、Mestre-Shanks算法、Schoof算法，其中Schoof算法效果不是很理想，一般用前两个。。。）
class Elliptic():
def __init__(self, a: int = 0, b: int = 0, char: int = 0, ):
'''
a is coefficient of EC
b is coefficient of EC
char is characteristic of finite field
'''
self.a = a
self.b = b
self.char = char
self.K = 0
delta = (4 * a ** 3 + 27 * b ** 2) % char
if delta == 0:
warnings.warn('delta is zero, may cause cusp form')
assert is_prime(char) == True, 'characteristic must be prime'
def __y2(self, x):
return (x ** 3 + self.a * x + self.b) % self.char
def get_y(self, x: int) -> False or int:
y2 = self.__y2(x)
if legendre(y2, self.char) == -1:
return False
else:
return quad_res(y2, self.char)
def point_x(self, x: int) -> (RationalPoint, RationalPoint) or False:
y = self.get_y(x)
if y:
return RationalPoint(x, y[0], self.char, self), RationalPoint(x, y[1], self.char, self)
else:
return False
def embed_plain(self, m: int) -> False or RationalPoint:
K = int(self.char / (m + 1))
t = 0
self.K = K
for i in range(K + 1):
x = m * self.K + i
y2 = self.__y2(x)
if legendre(y2, self.char) == 1:
t = 1
break
if t == 0:
return False
else:
y = self.get_y(x)
if type(y) is tuple:
if abs(m - y[0]) < abs(m - y[1]):
return RationalPoint(x, y[0], self.char, self)
else:
return RationalPoint(x, y[1], self.char, self)
else:
return RationalPoint(x, y, self.char, self)
def reco_plain(self, p: RationalPoint) -> False or int:
if self.K == 0:
return False
return p.x // self.K
def __str__(self):
s = 'y^2 = x^3 +' + str(self.a) + 'x+ ' + str(self.b) + ', charF=' + str(self.char)
return s
def __repr__(self):
s = 'y^2 = x^3 +' + str(self.a) + 'x+ ' + str(self.b) + ' in charF=' + str(self.char)
return s
def division_poly(self, n):
r = self.__gendivisionpoly(n)
return r
def __gendivisionpoly(self, n: int) -> list :
'''index of list = k+1'''
assert isinstance(n, int), 'the type of input is wrong'
if os.path.exists('division/')==False:
os.makedirs(r'division/')
root = 'division//'
binary_write(root+'f-1.pk',{0: -1})
binary_write(root+'f0.pk',{0: 0})
binary_write(root+'f1.pk',{0: 1})
binary_write(root+'f2.pk',{0: 2})
binary_write(root+'f3.pk',{0: -self.a**2, 1: 12 * self.b, 2: 6 * self.a, 4: 3})
binary_write(root+'f4.pk',(Polynomials({0: (-1)*self.a**3-8*self.b**2, 1: -4 *self.a * self.b,
2:-5*self.a**2,3:20*self.b, 4: 5*self.a,6:1},self.char)*4).ind_coe)
next_index = 5
k = n - 4
while True:
if k<=0:
break
if next_index%2==0:
t = (next_index) // 2
fn = Polynomials(binary_read(root+'f'+str(t)+'.pk'),self.char)
fa2 = Polynomials(binary_read(root+'f'+str(t+2)+'.pk'),self.char)
fa1 = Polynomials(binary_read(root+'f'+str(t+1)+'.pk'),self.char)
fm1 = Polynomials(binary_read(root+'f'+str(t-1)+'.pk'),self.char)
fm2 = Polynomials(binary_read(root+'f'+str(t-2)+'.pk'),self.char)
rr = multi_Nmul(fn,fa2,fm1,fm1)-multi_Nmul(fm2,fa1,fa1)
rr = rr%self.char
binary_write(root + 'f' + str(t * 2 ) + '.pk',rr.ind_coe)
else:
t = (next_index-1)//2
fn = Polynomials(binary_read(root+'f'+str(t)+'.pk'),self.char)
fa2 = Polynomials(binary_read(root+'f'+str(t+2)+'.pk'),self.char)
fa1 = Polynomials(binary_read(root+'f'+str(t+1)+'.pk'),self.char)
fm1 = Polynomials(binary_read(root+'f'+str(t-1)+'.pk'),self.char)
y = 4*Polynomials({3:1,1:self.a,0:self.b},self.char)
if next_index%2==0:
rr = multi_Nmul(y,y,fa2,fn,fn,fn)-multi_Nmul(fa1,fa1,fa1,fm1)
rr = rr % self.char
binary_write(root+'f'+str(t*2+1)+'.pk',rr.ind_coe)
else:
rr = multi_Nmul(fa2,fn,fn,fn)-multi_Nmul(fa1,fa1,fa1,fm1,y,y)
rr = rr % self.char
binary_write(root + 'f' + str(t * 2 + 1) + '.pk', rr.ind_coe)
next_index += 1
k -= 1
def __naive_counting(self):
t =0
p = self.char
for x in range(p):
c = self.point_x(x)
if c:
if c[0]==c[1]:
t+=1
else:
t += len(c)
return t+1
def cal_ord(self,P):
"""BSGS shanks counting"""
p = self.char
Q = (p+1)*P
m = int(pow(p,1/4)+1)
O = RationalPoint(0, 0, self.char, self)
for k in range(-m,m+1):
PP = Q + k * (2 * m * P)
t = 0
key = False
for item in PointList(P,m):
if PP==item:
t = t
key = True
break
elif -1*PP==item:
t = -t
key = True
break
t += 1
if key == True:
j = t
break
M = p+1+2*m*k-j
factor = factor_num(M)
t = 0
while True:
pi = list(factor.keys())[t]
cond = (M//pi)*P
if cond == O:
if M//pi==0:
break
M = M // pi
continue
t+=1
if t>= len(factor.keys()) :
break
return M
def cal_rationalpoint(self):
O = RationalPoint(0,0, self.char, self)
p = self.char
if p<=299:
return self.__naive_counting()
t = 1
lcm_M = []
lowerbound = p+1-2*math.sqrt(p)
upperbound = p+1+2*math.sqrt(p)
while True:
P = self.point_x(t)
if P:
P = P[0]
M = self.cal_ord(P)
lcm_M.append(M)
mlcm = multi_lcm(lcm_M)
if mlcm>=lowerbound and mlcm<=upperbound:
break
if t>=200:
break
t += 1
return mlcm
def schoof(self):
f_k = lambda k:Polynomials(binary_read('division//f'+str(k)+'.pk'),self.char)
factorlist = []
M = 1
t = 0
print('x')
while True:
if M>4*math.sqrt(self.char):
break
else:
p = prime_list[t]
M = M*p
factorlist.append(p)
t += 1
q = self.char
y = Polynomials({3:1,1:self.a,0:self.b},q)
t_list = []
g = y.gcd_poly(Polynomials({q:1,1:-1},q))
if g.leading_ind>0:
t_list.append((0,2))
else:
t_list.append((1,2))
print(factorlist)
for li in factorlist[1:]:
k = q%li
print(li)
xq = Polynomials({q**2:1,1:-1},q)
if k &1:
result = (sparse_multi(xq,f_k(k)*f_k(k))+multi_Nmul(f_k(k-1),f_k(k+1),y))
result.char = self.char
r1 = result.gcd_poly(f_k(li))
else:
result = (sparse_multi(xq, f_k(k)*f_k(k)*y) + multi_Nmul(f_k(k - 1), f_k(k + 1)))
result.char = self.char
r1 = result.gcd_poly(f_k(li))
if r1 != Polynomials({0:1}):
rk = legendre(q,li)
if rk!=-1:
pw,nw = quad_res(q,li)
xqq = Polynomials({q: 1, 1: -1},q)
if pw&1:
result = sparse_multi(xqq,multi_Nmul(f_k(pw),f_k(pw)))+multi_Nmul(f_k(pw-1),f_k(pw+1),y)
result.char = self.char
r2 = result.gcd_poly(f_k(li))
else:
result = sparse_multi(xqq,multi_Nmul( f_k(pw), f_k(pw),y)) + multi_Nmul(f_k(pw - 1), f_k(pw + 1))
result.char = self.char
r2 = result.gcd_poly(f_k(li))
if r2 != Polynomials({0:1}):
if pw&1:
result = 4*multi_Nmul(Npow_poly(y,(q-1)//2),f_k(pw),f_k(pw))-
multi_Nmul(f_k(pw-2),f_k(pw-1),f_k(pw-1))+multi_Nmul(f_k(pw-2),f_k(pw+1),f_k(pw+1))
result.char = self.char
r3 = result.gcd_poly(f_k(li))
else:
result = 4 * multi_Nmul(Npow_poly(y, (q + 3) // 2), f_k(pw), f_k(pw)) -
multi_Nmul(f_k(pw - 2), f_k(pw - 1), f_k(pw - 1)) + multi_Nmul(f_k(pw - 2), f_k(pw + 1), f_k(pw + 1))
result.char = self.char
r3 = result.gcd_poly(f_k(li))
if r3 == Polynomials({0:1}):
t = 2*nw %li
else:
t = 2*pw%li
else:
t = 0
else:
a3 = Polynomials({q**2:1,1:-1},q)
print('last')
if k&1:
a1 = multi_Nmul(y,f_k(k-1),f_k(k-1))
a2 = multi_Nmul(f_k(k),f_k(k))
a4 = multi_Nmul(a3,a2)-a1
a4 = multi_Nmul(a4,a4)
a5 = multi_Nmul(f_k(k+2),f_k(k-1),f_k(k-1))-multi_Nmul(f_k(k-2),f_k(k+1),f_k(k+1))
alpha2 = Npow_poly(multi_Nmul(Npow_poly(y,2),a5-4*Npow_poly(y,(q**2-1)//2)*Npow_poly(f_k(k),3)),2)
beta = 16*multi_Nmul(a2,a4,y)
else:
a1 = multi_Nmul(f_k(k-1),f_k(k-1))
a2 = multi_Nmul(y,f_k(k),f_k(k))
a4 = multi_Nmul(a3 , a2) - a1
a4 = multi_Nmul(a4,a4)
alpha2 = Npow_poly(
multi_Nmul(y, a5 - 4 * Npow_poly(y, (q ** 2 +3) // 2) * Npow_poly(f_k(k), 3)), 2)
beta = 16 * multi_Nmul(a2, a4, y)
rr1 = (multi_Nmul(a1-sparse_multi(Polynomials({q**2:1,q:1,1:1},q),a2),beta)+a2*alpha2)%f_k(li)
rr2 = 16*multi_Nmul(y,a2,beta)%f_k(li)
for item in range(1,li):
if item&1:
psi_tau = f_k(item)%f_k(li)
psi_tau = Npow_poly(psi_tau,q*2)%f_k(li)
psi_taupm = (y*(f_k(item+1)*f_k(item-1)))%f_k(li)
psi_taupm = Npow_poly(psi_taupm,q)%f_k(li)
else:
psi_tau = multi_Nmul(f_k(item),f_k(item),y)% f_k(li)
psi_tau = Npow_poly(psi_tau, q ) % f_k(li)
psi_taupm = (f_k(item + 1) * f_k(item - 1)) % f_k(li)
psi_taupm = Npow_poly(psi_taupm, q) % f_k(li)
if (rr1*psi_tau+rr2*psi_taupm)%f_k(li)==Polynomials({0:0},q):
t = item%li
break
if t==li:
t -= 1
break
t_list.append((t,li))
res,mod =[],[]
for i in t_list:
res.append(i[0])
mod.append(i[1])
result = CRT(res,mod)
M = multi_lcm(mod)
while True:
if result>q+1-np.sqrt(q):
break
result = result + M
return result
再补两个函数，这两个是用于Schoof算法中，降低内存开销的策略。
def binary_write(filename,ob):
with open(filename,'wb') as f:
pk.dump(ob,f)
def binary_read(filename):
with open(filename,'rb') as f:
result = pk.load(f)
return result
numpy复数类的重载
因为原始的库在计算大数据时会出现类型错误（也有可能我这个包有点过时了。。。。），所以我对此进行了一定的修改，并允许“代数整数”的定义。
class Complex(np.complex):
def __add__(self, other):
r = super().__add__(other)
return Complex(r.real, r.imag)
def __radd__(self, other):
return self+other
def __mul__(self, other):
r = super().__mul__(other)
return Complex(r.real, r.imag)
def __sub__(self, other):
return self+(-1)*other
def __pow__(self, power, modulo=None):
r = super().__pow__(power)
return Complex(r.real, r.imag)
def __rtruediv__(self, other):
return Complex(super().__rtruediv__(other).real, super().__rtruediv__(other).imag)
def __rdiv__(self, other):
r = super().__truediv__(other)
return Complex(r.real, r.imag)
def __mod__(self, other):
assert other!=0,'divided by zero'
return Complex(self.real%other,self.imag%other)
def __truediv__(self, other):
r = super().__truediv__(other)
return Complex(r.real, r.imag)
def __round__(self, n=None):
return Complex(round(self.real,n),round(self.imag,n))
有理点类
实现了有理点中的 加、减、数乘
其中，加法是主要编写的，数乘的算法借用了快速幂的思想。
class RationalPoint():
'''for elliptic curve in finite field'''
def __init__(self, x, y, char, e):
self.x = x%char
self.y = y%char
self.char = char
self.E = e
def __sub__(self, other):
p = RationalPoint(other.x, self.char - other.y, other.char, other.E)
q = RationalPoint(self.x, self.y, self.char, self.E)
return q + p
def __mul__(self, other):
if type(other) == int:
r = RationalPoint(self.x, self.y, self.char, self.E)
if other<0:
r = RationalPoint(self.x, - self.y, self.char, self.E)
other = -other
if other == 0:
return RationalPoint(0, 0, self.char, self.E)
elif other == 1:
return r
ans = RationalPoint(0, 0, self.char, self.E)
while other!=0:
if other&1:
ans = ans+r
r = r+r
other>>=1
return ans
def __neg__(self):
return -1*self
def __rmul__(self, other):
r = self*other
return r
def __add__(self, other):
p = RationalPoint(other.x, other.y, other.char, other.E)
x_1, y_1 = self.x, self.y
x_2, y_2 = p.x, p.y
char = self.char
if x_2 == 0 and y_2 == 0:
return self
elif x_1 == 0 and y_1 == 0:
return other
if x_1 == x_2 and y_1 != y_2:
return RationalPoint(0, 0, self.char, self.E)
elif x_1 != x_2:
k = (y_2 - y_1) * inv_fin(x_2 - x_1, char)
b = y_1 - k * x_1
x_3 = k ** 2 - x_1 - x_2
x_3 = x_3 % char
y_3 = -(k * x_3 + b)
return RationalPoint(x_3 % char, y_3 % char, self.char, self.E)
else:
if y_2==0:
return RationalPoint(0, 0, self.char, self.E)
x_3 = (((3 * x_1 ** 2 + self.E.a) * inv_fin(2 * y_1, char)) ** 2 - 2 * x_1)
y_3 = -y_1 + ((3 * x_1 ** 2 + self.E.a) * inv_fin(2 * y_1, char)) * (x_1 - x_3)
return RationalPoint(x_3 % char, y_3 % char, self.char, self.E)
def __str__(self):
s = '(' + str(self.x) + ',' + str(self.y) + ')'
return s
def __repr__(self):
s = '(' + str(self.x) + ',' + str(self.y) + ')'
return s
def __eq__(self, other):
if self.x == other.x and self.y == other.y:
return True
else:
return False
最大公因数
def gcd(a: int, b: int) -> int:
assert isinstance(a, int) and isinstance(b, int), 'the type of input is wrong'
while True:
if a % b == 0:
return b
else:
a, b = b, a % b最小公倍数
def lcm(a: int, b: int) -> int:
assert isinstance(a, int) and isinstance(b, int), 'the type of input is wrong'
return a * b // gcd(a, b)数组的最大公因数
def multi_lcm(a: list) -> int:
result = [i for i in a]
for i in range(len(result)-1):
lcm_tmp = lcm(result[i],result[i+1])
result[i+1] = lcm_tmp
return result[-1]数组的最小公倍数
def multi_gcd(a: list) -> int:
assert isinstance(a, list), 'the type of input is wrong'
r = sorted(a, reverse=True)
r.append(r[0])
while True:
if len(set(r[:-1])) == 1:
break
for i in range(len(r) - 1):
if r[i] % r[i + 1] == 0:
r[i] = abs(r[i + 1])
else:
r[i] = abs(r[i] % r[i + 1])
r[-1] = r[0]
return r[0]
穷举法生成特定点数的椭圆曲线
属于暴力穷举法
def finEC(N,p):
"""naive method"""
x = 1
t = p+1-N
k = 0
for item in range(p):
if legendre(item,p)==-1 :
k = item
break
for i in range(p):
if 4*i**3+27*i**2%p!=0:
Ea = Elliptic(i,-i,p)
P = Ea.point_x(x)
if (p+1-t)*P[0]==0*P[0]:
return Ea
elif (p+1+t)*P[0]==0*P[0]:
Eb = Elliptic((k**2)*i,-(k**3)*i,p)
return Eb
return False类数算法
其实类数算法有很多，这里的算法是为了类域多项式的估计而写的。
def hei_D(D:int )->int :
assert D%4 in [0,1],'not fundamental discriminant'
h = 1
b = D%2
B = int(math.sqrt(abs(D)/3))
while True:
q = (b**2-D)//4
a = b
if a <1 :
a = 1
while True:
if q%a==0:
if a==b or a**2==q or b==0:
h+=1
else:
h+=2
a += 1
if a**2<=q:
continue
while True:
b += 2
if b<=B:
break
else:
return h-1
breakHilbert 类域多项式的生成
def getHCP(D,prime):
#Hilbert class polynomial normal generation
assert D % 4 in [0, 1], 'not fundamental discriminant'
p = Polynomials({0:1})
b = D%2
B = int(math.sqrt(abs(D)/3))
while b<=B:
t = (b ** 2 - D) // 4
a = max(b, 1)
while a**2<=t:
if t%a==0:
j = j_invar((-b+np.power(Complex(D,0),0.5))/(2*a))
j = Complex(j.real,j.imag)
if a==b or a**2==t or b==0:
p = sparse_multi(p,Polynomials({1:1,0:-j}))
else:
p =sparse_multi(p,Polynomials({2:1,1:-2*j.real,0:j.real**2+j.imag**2}))
a += 1
b += 2
if b>B:
ind_coe = {i:round(Complex(p.ind_coe[i].real,p.ind_coe[i].imag),0) for i in p.ind_coe}
for item in ind_coe:
if ind_coe[item].imag==0:
ind_coe[item] = ind_coe[item].real
ind_coe[item] = int(ind_coe[item])
return Polynomials(ind_coe,prime)%prime
j不变量的数值计算
def j_invar(tau):
"""calculate the j-invariant via tau"""
M = 19
m = np.arange(1, M+1)
left_part = np.ones(m.shape, dtype=np.float)
left_part[m%2 == 1] = -1
exponent = 0.5*m*(3*m-1)
exponent2 = 0.5*m*(3*m+1)
q = lambda t: np.e ** (2 * np.pi * Complex(0, 1) * t)
q = np.vectorize(q)
right_part1 = q(tau)**exponent
right_part2 = q(tau)**exponent2
right_part12 = q(2*tau)**exponent
right_part22 = q(2 * tau) ** exponent2
eta = q(tau/24)*(1+np.sum(left_part*(right_part1+right_part2)))
eta2 = q(tau/12)*(1+np.sum(left_part*(right_part12+right_part22)))
ft = np.float_power(eta2 / eta , 24)
j = np.float_power((256 * ft + 1),3) / ft
return j曲线生成（CM法）
def gen_EC( N,p):
t = p+1-N
D = 4*p-t**2
D = squarefreefactor(D)
H_D = getHCP(-D,p)
j = 1
while True:
if H_D.to_lambda()(j)%p!=0:
j+=1
else:
break
if j >=p:
print('no solve')
a = 1
k = j*inv_fin(1728-j,p)
e = Elliptic(3*k,2*k,p)
i = 2
while True:
P = e.point_x(i)
if P:
P = P[0]
break
i += 1
if N*P!=0*P:
c = 0
for item in range(p):
if legendre(item, p) == -1:
c = item
break
e = Elliptic((c ** 2) *3*k%p, (c ** 3) * 2*k%p, p)
return e基于CRT算法的Hilbert类域多项式计算并生成曲线
和上一个类比，这个算法主要是用于规模足够大的有限域特征的
不过这个算法是基于有理点个数计算的。。。目前还没有实现比较快速的算法。。
def modCRT(n: int ,Sm: list,Sa: list ,epsi: float) -> int :
M = 1
for i in Sm:
M *= i
Mi = []
ai = []
z = 0
for i in range(len(Sm)):
Mi_tmp = M//Sm[i]
ai_tmp = inv_fin(Mi_tmp,Sm[i])%Sm[i]
z += ai_tmp*Mi_tmp*Sa[i]
Mi.append(Mi_tmp)
ai.append(ai_tmp)
r = int(z/M+0.5)%n
M = M%n
z = z %n
return (z-r*M)%n
def HCP_CRT(D,p):
assert D % 4 in [0, 1], 'not fundamental discriminant'
B = int(math.sqrt(abs(D) / 3))
b = 0
a_lst = []
while b <= B:
if (b**2-D)%4!=0:
b+=1
continue
t = (b ** 2 - D) // 4
a = max(b, 1)
while a ** 2 <= t:
if t %a!=0:
a += 1
continue
c = t//a
print(a,b,c,'s')
if a==b or a**2==(b**2-D)//4 or b==0:
a_lst.append(1 / a)
else:
a_lst.append(1/a)
a_lst.append(-1/a)
a+=1
b += 1
if b > B:
break
h = hei_D(D)
Bound = comb(h,int(h/2))*np.exp(np.pi*np.sqrt(abs(D))*np.sum(a_lst))
S = []
H = []
M = 1
t = 1
while M<=Bound:
qf = t**2+abs(D)
if qf%4!=0:
t += 1
continue
else:
q = qf//4
if is_prime(q):
S.append(q)
M = M*q
t += 1
for q in S:
Sq = []
t = 0
for j in range(q):
if j==0 or j==1728:
continue
k = j*inv_fin(1728-j,q)
E = Elliptic(3*k%q,2*k%q,q)
points = E.cal_rationalpoint()
tt = abs(np.sqrt(4 * q - abs(D)))
if points==q+1+tt or points== q+1-tt:
Sq.append(j)
print(E,j)
t+=1
if t>=h:
break
Hq = Polynomials({0:1})
for item in Sq:
Hq = sparse_multi(Hq,Polynomials({1:1,0:-item}))%q
H.append(Hq)
print('abcd')
c_k = []
for i in range(0,h):
Sa = [j[i] for j in H ]
c_k.append(modCRT(p,S,Sa,0.001))
r = {i: c_k[i] for i in range(len(c_k))}
r.update({h: 1})
result = Polynomials(r)
return result
有限域上多项式计算
参考如下文章吧
Alepha E：杂谈 快速数论变换 以及代数数实现上的思考zhuanlan.zhihu.com就这些算法了
扫一扫
963
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
