NTRU是目前据我所知的最快的抗量子加密算法(当然NTRU也有也有签名算法)。

关于算法的Introduction我就不多说了,如果你之前了解过NTRU,那么我写的不一定会比你所了解的多;如果没了解过,可以先去谷歌或维基了解一下,也可以参考其官网上的文章。

NTRU到目前为止已经有很多版本,我参考的是[HPS14] 7.10 节的教科书版本,顺带一题,作者也是提出NTRU的三人(https://ntru.org/f/hps96.pdf)

参考:

[HPS14] Hoffstein J, Pipher J, Silverman J H, et al. An introduction to mathematical cryptography[M]. New York: springer, 2008.

PS:第一版好像是08年的,但我参考的是14年的第二版。

[HPS17] Hoffstein J, Pipher J, Schanck J M, et al. Choosing parameters for NTRUEncrypt[C]//Topics in Cryptology–CT-RSA 2017: The Cryptographers’ Track at the RSA Conference 2017, San Francisco, CA, USA, February 14–17, 2017, Proceedings. Springer International Publishing, 2017: 3-18.

PS:这篇主要讲的是参数设置。

要看懂下面的内容你可能需要:

  • 数论知识(特别群论、环论);
  • 线性代数知识,多项式运算等;
  • 基础的密码学知识(起码你得知道加密算法是干什么的吧);
  • Sagemath代码的读写经验;
  • 足够的耐心。

理论知识·

多项式环和参数选取·

首先[HPS14]的7.10.1节一开始就扔出了三个卷积多项式环(Convolution Polynomial Rings),RRRpR_pRqR_q,理解这三个环其实是理解NTRU的关键,接下来细说。

首先看RR,分子(姑且先叫分子吧)的Z[x]\mathbb{Z}[x]指的是整数域上的多项式,以符号xx作为多项式的未知数,实际上也就是小学/初中开始接触的多项,比如x2+x+1x^2+x+12x3x-2x^3-x等。

然后,横线其实不是指除法,比起除法其实更像一个取余,或者叫商环(Quotient Ring,可以参考[HPS14]的2.10),大概意思是,“分子”上的多项式模上xN1x^N-1后就会落在RR上,其中N是NTRU的参数。其实跟RSA在Zn\mathbb{Z}_n^*中运算所以要模nn类似,可以做个类比方便理解。

举个栗子,假设取N=2N=2,那就是R=Z[x]x21R = \frac{\mathbb{Z}[x]}{x^2-1},现在尝试把Z[x]\mathbb{Z}[x]上的x2+x+1x^2+x+1转换到RR上,由于“模数”是x21x^2-1,所以有

x210 (in R)\begin{aligned} x^2-1 &\equiv 0 \ (\text{in } R) \end{aligned}

把这个关系利用到x2+x+1x^2+x+1上可以得到:

x2+x+1x2+x+10x2+x+1(x21)x+2 (in R)\begin{aligned} x^2+x+1 &\equiv x^2+x+1 - 0 \\ &\equiv x^2+x+1 - (x^2-1) \\ &\equiv x + 2 \ (\text{in } R) \end{aligned}

x+2x+2就是最终落在RR上的结果。

另一个栗子,2x3x-2x^3-x,这里偷一下懒,对其作分解后可以得到:

2x3x=2(x21)x3x-2x^3-x = -2(x^2-1)·x - 3x

所以有:

2x3x2(x21)x3x203x3x (in R)\begin{aligned} -2x^3-x &\equiv -2(x^2-1)·x - 3x \\ &\equiv -2·0 - 3x \\ &\equiv -3x \ (\text{in } R) \end{aligned}

一般来说RR中的多项式最高项不会大于等于NN次,因为xN1 (in R)x^N \equiv 1 \ (\text{in } R),所以每出现xNx^N都可以用11代替,可以利用这个方法检验最终结果的正确性(虽然我都是直接用Sage算的,逃)。

RpR_p其实类似,只不过其“分子”(Z/pZ)[x](\mathbb{Z}/p\mathbb{Z})[x]是系数模pp的多项式,也就是多项式的系数需要进行模pp操作,使其落在Zp\mathbb{Z}_p中(Z/pZ\mathbb{Z}/p\mathbb{Z}其实就是Zp\mathbb{Z}_p的另一种写法)。

举个栗子,设p=3p=3,下面尝试把2x3x-2x^3-x转换到(Z/pZ)[x](\mathbb{Z}/p\mathbb{Z})[x]中:

2x3x(2(modp))x3+(1(modp))x(2(mod3))x3+(1(mod3))xx3+2x (in (Z/pZ)[x])\begin{aligned} -2x^3-x &\equiv (-2 \pmod p)x^3 + (-1 \pmod p)x \\ &\equiv (-2 \pmod 3)x^3 + (-1 \pmod 3)x \\ &\equiv x^3 + 2x \ (\text{in } (\mathbb{Z}/p\mathbb{Z})[x]) \end{aligned}

其实只是系数模pp

然后RpR_p中的取余操作和RR的类似,不再累赘。RqR_qRpR_p类似,只是模数不同。

以上的NNppqq是NTRU的公开参数,在[HPS14]中的参数要求是,NN为素数,且

gcd(N,q)=gcd(p,q)=1\rm{gcd}(N, q) = \rm{gcd}(p, q) = 1

以下更详细的参数设置建议摘抄自[HPS17],大概意思是:pp一般取p=3p=3qq为了加快运算一般取22的幂次方(盲猜是二进制有利于计算机优化),NN同样取素数,不过某些特殊版本的基于理想格的NTRU会取形如N=2nN=2^n,然后模数取xN+1x^N+1,这里超纲了,略。

三元多项式和密钥生成·

接下来会说到一个叫T\mathcal{T}的函数:

T\mathcal{T}的输入为两个整数d1d_1d2d_2,输出为一个三元多项式(Ternary Polynomial),三元即(1,0,1)(-1, 0, 1),大概意思是,T\mathcal{T}函数会采样一个符合条件:有d1d_1个项系数为11、有d2d_2个项系数为1-1,其余项系数为00(即没有)的多项式,且这个多项式落在RR上。

接下来利用T\mathcal{T}采样密钥,首先是私钥,我就直接上图不再码一遍了:

其中的dd也是公开参数(上面也有说了),f(x)f(x)g(x)g(x)是两个主要的私钥(甚至可以把g(x)g(x)看作一个随机数),落在RR中(因为T\mathcal{T}的返回值落在RR中,先记住,后面会用到),Fq(x)F_q(x)Fq(x)F_q(x)分别是f(x)f(x)RqR_qRpR_p上的逆元,这里写出来的意思是,把Fq(x)F_q(x)Fq(x)F_q(x)存储起来可以加快运算速度,因为计算多项式环的逆元需要一段时间(虽然也不是很久)。

多嘴说一句,f(x)f(x)取自T(d+1,d)\mathcal{T}(d+1, d)是因为需要f(x)f(x)可逆,书中有提到T(d,d)\mathcal{T}(d, d)的不可逆,也可以自己推一遍。

这里隐藏的一个信息是,RR上的多项式可以转换到RqR_qRpR_p,系数模qq或者模pp而已,挺显然的,略。

然后计算公钥

换一种写法其实也就是:

h(x)f(x)1g(x) in Rqh(x) \equiv f(x)^{-1} · g(x) \ \text{in } R_q

Center-lift和加解密·

首先说Center-lift,上面说了,RR上的多项式可以转换到RqR_qRpR_p中,而Center-lift则是把RqR_qRpR_p中的多项式转换回RR中,而且使得多项式系数的绝对值最小(相比于其他的转换方式)。

Center-lift的定义可以看[HSP14]的7.9节,主要是针对多项式系数的操作,可以和补码类比。

上面假设是RqR_q多项式的Center-lift,即RqR_q转换到RR。上图的a(x)a(x)RqR_q上的一个多项式,a(x)a'(x)为其Center-lift到RR后的结果,aia_ia(x)a(x)ii次项系数,aia_i'a(x)a'(x)ii次项系数。

Center-lift后的结果是,所有的系数aia_i'会落在区间(q2,q2](-\frac{q}{2}, \frac{q}{2}],而且还能保证aiai(modq)a_i' \equiv a_i \pmod q

Center-lift的方法也不难,和转补码类似,在q2\frac{q}{2}砍一半,小于等于q2\frac{q}{2}的不动,大于q2\frac{q}{2}的平移到(q2,0)(-\frac{q}{2}, 0),也就是:

ai={ai,if 0aiq2aiq,if q2<ai<qa_i' = \begin{cases} a_i, & \text {if $0 \le a_i \le \frac{q}{2}$} \\ a_i - q, & \text {if $\frac{q}{2} < a_i < q $} \end{cases}

检验一下,上面的划分可以覆盖整个Zq\mathbb{Z}_q,当0aiq20 \le a_i \le \frac{q}{2}时显然aiai(modq)a_i' \equiv a_i \pmod q,而且落在(q2,q2](-\frac{q}{2}, \frac{q}{2}]q2<ai<q\frac{q}{2} < a_i < q时,aiaiqai0ai(modq)a_i' \equiv a_i - q \equiv a_i - 0 \equiv a_i \pmod q,由于aia_i落在(q2,q)(\frac{q}{2}, q),所以aiaiqa_i' \equiv a_i - q落在(q2q,qq)(\frac{q}{2}-q, q-q),也即(q2,0)(-\frac{q}{2}, 0),自然落在(q2,q2](-\frac{q}{2}, \frac{q}{2}],检验通过。

然后接着看加密,明文空间是RpR_p做Center-lift后的空间,也即在RR上,系数取值落在(p2,p2](-\frac{p}{2}, \frac{p}{2}]的空间,加密使用的随机数r(x)r(x)T(d,d)\mathcal{T}(d, d)中采样,然后加密操作就是下图的公式,密文是e(x)e(x)

在继续解密操作前,先分析一下密文的结构。

首先如果要解密的话,其实只用把ph(x)r(x)p·h(x)·r(x)这一项去除,先剧透一下,由于p0 in Rpp \equiv 0 \ \text{in } R_p,所以直观上把e(x)e(x)放入RpR_p中就可以把ph(x)r(x)p·h(x)·r(x)去除,而由于明文空间是RpR_p做Center-lift后的空间,所以把m(x) in Rqm(x) \ \text{in }R_q放入RpR_p中并不会改变m(x)m(x)的结构,即可以恢复明文。

剩下的问题是我其实并不能保证[e(x)(modq)]e(x)(modp)[e(x) \pmod q] \equiv e(x) \pmod p,除非q=kpq = kp(这种情况在前几天讲OU98时有讲过),但因为gcd(p,q)=1\rm{gcd}(p, q) = 1,所以这种情况并不存在。

这样的栗子可以举很多,可以直接看整数的栗子(反正也是模在多项式系数上),设e=65537e=65537q=512q=512p=3p=3,计算可得

(e(modq))(modp)=(65537(mod512))(mod3)=1(mod3)=1e(modp)=65537(mod3)=2\begin{aligned} (e \pmod q) \pmod p = (65537 \pmod {512}) \pmod 3 = 1 \pmod 3 &= 1 \\ e \pmod p = 65537 \pmod 3 &= 2 \end{aligned}

显然不等,也即RqR_q不能直接转换到RpR_p上。

但是前面说过,RR可以转换到RpR_p上,所以可以考虑先把e(x)e(x)转换(Center-lift)到RR上,然后再转到RpR_p上。

但新的问题是,直接把RqR_q上的e(x)e(x)转到RR上,其结果不一定是原来RR上的e(x)e(x),举个整数域上的栗子即65537(mod512)=16553765537 \pmod {512} = 1 \ne 65537,除了一种情况,即把e(x)e(x)放在RR上运算后的系数本来就落在(q2,q2](-\frac{q}{2}, \frac{q}{2}],也即RqR_q上的元素做Center-lift后的区间。

在讲如何把e(x)e(x)系数转化到(q2,q2](-\frac{q}{2}, \frac{q}{2}]前,可能先需要了解一下多项式的乘法。

多项式卷积乘法·

在继续讲解密思路前,先“简要”介绍一下多项式卷积。

这篇我曾经写过,多项式的卷积乘法可以看成一个矩阵操作,这里多嘴推导一下这个矩阵是怎么来的,如果不关心推导则可以跳过

下面假设f(x)f(x)g(x)g(x)RR中相乘(假设而已),先把多项式表述为:

f(x)i=0N1fixig(x)i=0N1gixif(x) \equiv \sum^{N-1}_{i=0} f_i x^i \\ g(x) \equiv \sum^{N-1}_{i=0} g_i x^i

那么有(参考多项式乘法,用分配律推也行):

f(x)g(x)i=0N1(fixij=0N1gjxj)i=0N1j=0N1fixigjxji=0N1j=0N1figjxi+j (in R)\begin{aligned} f(x) · g(x) &\equiv \sum^{N-1}_{i=0} (f_i x^i \sum^{N-1}_{j=0} g_j x^j) \\ &\equiv \sum^{N-1}_{i=0} \sum^{N-1}_{j=0} f_i x^i g_j x^j \\ &\equiv \sum^{N-1}_{i=0} \sum^{N-1}_{j=0} f_i g_j x^{i+j} \ \text{(in $R$)} \end{aligned}

由于f(x)g(x)f(x) · g(x)的结果肯定也是个多项式(封闭性),所以先把上面结构调整到像一个多项式,即合并次数相同的项(注意f(x)f(x)g(x)g(x)最高项只能是xN1x^{N-1},所以不考虑模的时候乘起来最高项只能是xN1xN1=x2(N1)x^{N-1} · x^{N-1} = x^{2(N - 1)},然后小于NN次的项和大于等于NN次的项分开处理):

f(x)g(x)i=0N1j=0N1figjxi+jk=0N1((i=0kfigki)xk)+k=N2(N1)((i=k(N1)N1figki)xk) (in R)\begin{aligned} f(x) · g(x) &\equiv \sum^{N-1}_{i=0} \sum^{N-1}_{j=0} f_i g_j x^{i+j} \\ \\ &\equiv \sum^{N-1}_{k=0} ((\sum^{k}_{i=0} f_{i} g_{k-i}) · x^k) + \sum^{2(N-1)}_{k=N} ((\sum^{N-1}_{i=k-(N-1)} f_{i} g_{k-i}) · x^k) \ \text{(in $R$)} \end{aligned}

最开始讲环的时候说到,RR是模xN1x^N-1的,所以xN1x0 (in R)x^N \equiv 1 \equiv x^0 \ \text{(in }R \text{)},所以有

xN+ixi (in R)x^{N + i} \equiv x^i \ \text{(in $R$)}

即未知数xx的指数部分需要模NN,所以可以继续化简:

f(x)g(x)k=0N1((i=0kfigki)xk)+k=N2(N1)((i=k(N1)N1figki)xkN)k=0N1((i=0kfigki)xk)+k=0N2((i=(k+N)(N1)N1fig(k+N)i)xk)k=0N1((i=0kfigki)xk)+k=0N1((i=k+1N1figk+Ni)xk)k=0N1(i=0kfigki+i=k+1N1figk+Ni)xkk=0N1(i=0N1figki(modN))xk (in R)\begin{aligned} f(x) · g(x) &\equiv \sum^{N-1}_{k=0} ((\sum^{k}_{i=0} f_{i} g_{k-i}) · x^k) + \sum^{2(N-1)}_{k=N} ((\sum^{N-1}_{i=k-(N-1)} f_{i} g_{k-i}) · x^{k-N}) \\ &\equiv \sum^{N-1}_{k=0} ((\sum^{k}_{i=0} f_{i} g_{k-i}) · x^k) + \sum^{N-2}_{k=0} ((\sum^{N-1}_{i=(k+N)-(N-1)} f_{i} g_{(k+N)-i}) · x^k) \\ &\equiv \sum^{N-1}_{k=0} ((\sum^{k}_{i=0} f_{i} g_{k-i}) · x^k) + \sum^{N-1}_{k=0} ((\sum^{N-1}_{i=k+1} f_{i} g_{k+N-i}) · x^k) \\ &\equiv \sum^{N-1}_{k=0} (\sum^{k}_{i=0} f_{i} g_{k-i} + \sum^{N-1}_{i=k+1} f_{i} g_{k+N-i}) · x^k \\ &\equiv \sum^{N-1}_{k=0} (\sum^{N-1}_{i=0} f_{i} g_{k-i \pmod N}) · x^k \ \text{(in $R$)} \end{aligned}

其中的i=0N1figki(modN)\sum^{N-1}_{i=0} f_{i} g_{k-i \pmod N}即为f(x)g(x)f(x) · g(x)kk次项系数。

粗略检验一下正确性,设f(x)=x2+2x+3f(x)=x^2+2x + 3g(x)=3x2+2x+1g(x) = 3x^2 + 2x + 1N=3N=3,则:

f(x)g(x)(x2+2x+3)(3x2+2x+1)3x4+8x3+14x2+8x+33x+8+14x2+8x+314x2+11x+11 (in R)\begin{aligned} f(x) · g(x) &\equiv (x^2+2x + 3) · (3x^2 + 2x + 1) \\ &\equiv 3 x^4 + 8 x^3 + 14 x^2 + 8 x + 3 \\ &\equiv 3x + 8 + 14x^2 + 8x + 3 \\ &\equiv 14x^2 + 11x + 11 \ \text{(in $R$)} \end{aligned}

k=0N1(i=0N1figki(modN))xkk=02(i=02figki(mod3))xk(i=02fig0i(mod3))x0+(i=02fig1i(mod3))x1+(i=02fig2i(mod3))x2(f0g0+f1g2+f2g1)x0+(f0g1+f1g0+f2g2)x1+(f0g2+f1g1+f2g0)x2(31+23+12)x0+(32+21+13)x1+(33+22+11)x214x2+11x1+11x014x2+11x+11 (in R)\begin{aligned} &\sum^{N-1}_{k=0} (\sum^{N-1}_{i=0} f_{i} g_{k-i \pmod N}) · x^k \\ \equiv &\sum^{2}_{k=0} (\sum^{2}_{i=0} f_{i} g_{k-i \pmod 3}) · x^k \\ \equiv &(\sum^{2}_{i=0} f_{i} g_{0-i \pmod 3}) · x^0 + (\sum^{2}_{i=0} f_{i} g_{1-i \pmod 3}) · x^1 + (\sum^{2}_{i=0} f_{i} g_{2-i \pmod 3}) · x^2 \\ \equiv &(f_0 g_0 + f_1 g_2 + f_2 g_1) · x^0 + (f_0 g_1 + f_1 g_0 + f_2 g_2) · x^1 + (f_0 g_2 + f_1 g_1 + f_2 g_0) · x^2 \\ \equiv &(3·1 + 2·3 + 1·2) · x^0 + (3·2 + 2·1 + 1·3) · x^1 + (3·3 + 2·2 + 1·1) · x^2 \\ \equiv &14x^2 + 11x^1 + 11x^0 \\ \equiv &14x^2 + 11x + 11 \ \text{(in $R$)} \end{aligned}

检验通过。

以下用矩阵的形式表达k=0N1(i=0N1figki(modN))xk\sum^{N-1}_{k=0} (\sum^{N-1}_{i=0} f_{i} g_{k-i \pmod N}) · x^k,姑且记f(x)g(x)f(x)·g(x)ii次项系数为λi\lambda_i,即为

(f0,f1,...,fN1)(g0g1g2gN1gN1g0g1gN2gN2gN1g0gN3g1g2g3g0)=(λ0,λ1,...,λN1)  (in R)(f_0, f_1, ..., f_{N-1}) · \begin{pmatrix} g_0 & g_1 & g_2 & \cdots & g_{N-1} \\ g_{N-1} & g_0 & g_1 & \cdots & g_{N-2} \\ g_{N-2} & g_{N-1} & g_0 & \cdots & g_{N-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ g_1 & g_2 & g_3 & \cdots & g_0 \\ \end{pmatrix} \\ = (\lambda_0, \lambda_1, ..., \lambda_{N-1})\ \ \text{(in $R$)}

解密与参数关系计算·

首先直接贴解密过程,

公式推导我也懒得再码一遍了,直接截图:

现在的目标是要把e(x)e(x)的系数转化到(q2,q2](-\frac{q}{2}, \frac{q}{2}]区间,而观察e(x)e(x)的结构可发现,要达成这个目标最大的阻碍是h(x)h(x),因为h(x)f(x)1g(x) in Rqh(x) \equiv f(x)^{-1} · g(x) \ \text{in } R_q,其中的f(x)1f(x)^{-1}取逆操作会把系数“打乱”,从而产生很大(qq规模)的系数,所以需要想方法消去。

消去方法如上面a(x)a(x)的推导,乘上一个f(x)f(x),把h(x)h(x)转换成g(x)g(x)g(x)T(d,d)g(x) \in \mathcal{T}(d, d),其系数只会是小系数1-10011

虽然消去h(x)h(x),但其实还不能保证a(x)a(x)的绝对值系数严格落在(q2,q2](-\frac{q}{2}, \frac{q}{2}],所以接下来需要先计算a(x)a(x)的最大和最小系数,计算方法稍微引用下原文:

可以把a(x)a(x)分成两部分,即“+”分割的pg(x)r(x)p · g(x) · r(x)f(x)m(x)f(x) · m(x)

首先看pg(x)r(x)p · g(x) · r(x)的最大和最小系数,由于pp是常数,所以可以先看g(x)r(x)g(x) · r(x),翻看上面内容,g(x)g(x)r(x)r(x)都采样自T(d,d)\mathcal{T}(d, d),即只有dd个系数为11的项和dd个系数为1-1的项,其余项系数为00

上一节已经推出了,i=0N1girki(modN)\sum^{N-1}_{i=0} g_{i} r_{k-i \pmod N}g(x)r(x)g(x) · r(x)kk次项系数,由于girki(modN)g_{i} r_{k-i \pmod N}只会在gig_{i}rki(modN)r_{k-i \pmod N}都取11或者都取1-1时可以取得最大值11,所以在拥有最多的girki(modN)g_{i} r_{k-i \pmod N}出现最大值11的情况下g(x)r(x)g(x) · r(x)kk次项系数出现最大值,这个情况即是:值为11gig_{i}刚好匹配上值为11rki(modN)r_{k-i \pmod N}、值为1-1gig_{i}刚好匹配上值为1-1rki(modN)r_{k-i \pmod N}、且值为00gig_{i}刚好匹配上值为00rki(modN)r_{k-i \pmod N}的情况,产生的最大值是2d2d

再把常数pp考虑进去,即pg(x)r(x)p · g(x) · r(x)的最大值是p2dp·2d

类似地,当111-1匹配上的时候,girki(modN)g_{i} r_{k-i \pmod N}出现最小值1-1,即得pg(x)r(x)p · g(x) · r(x)的最小值是p2d-p·2d

再来看f(x)m(x)f(x) · m(x)得最大和最小系数,f(x)f(x)采样自T(d+1,d)\mathcal{T}(d+1, d),即有d+1d+1个项系数为11、有dd个项系数为1-1、其余项为00m(x)m(x)的系数落在(p2,p2](-\frac{p}{2}, \frac{p}{2}],即RpR_p做Center-lift后的区间。

类似地,i=0N1fimki(modN)\sum^{N-1}_{i=0} f_{i} m_{k-i \pmod N}f(x)m(x)f(x) · m(x)kk次项系数,fimki(modN)f_{i} m_{k-i \pmod N}只会在fif_{i}11mki(modN)m_{k-i \pmod N}p2\frac{p}{2}时和在fif_{i}1-1mki(modN)m_{k-i \pmod N}p2-\frac{p}{2}(虽然取不到,但可以用来计算界)时取得最大值p2\frac{p}{2},这些情况总共可以出现(d+1)+d(d+1)+d次(f(x)f(x)11的个数加上1-1的个数),所以f(x)m(x)f(x) · m(x)的系数最大值是(2d+1)p2(2d+1) · \frac{p}{2}

再类似地,在fif_{i}11mki(modN)m_{k-i \pmod N}p2-\frac{p}{2}时和在fif_{i}1-1mki(modN)m_{k-i \pmod N}p2\frac{p}{2}时,fimki(modN)f_{i} m_{k-i \pmod N}取得最小值p2-\frac{p}{2},总可出现(2d+1)(2d+1)次,所以所以f(x)m(x)f(x) · m(x)的系数最小值是(2d+1)p2-(2d+1) · \frac{p}{2}

最后a(x)a(x)的最大系数就是pg(x)r(x)p · g(x) · r(x)f(x)m(x)f(x) · m(x)两部分的最大系数之和,即

p2d+(2d+1)p2=3dp+p2=(3d+12)pp·2d + (2d+1) · \frac{p}{2} = 3dp + \frac{p}{2} = (3d + \frac{1}{2}) · p

同理,最小系数是(3d+12)p-(3d + \frac{1}{2}) · p。(希望到这里还没看晕,逃)

接下来,要保证a(x)a(x)的系数落在(q2,q2](-\frac{q}{2}, \frac{q}{2}],只需要保证其最大系数严格小于q2\frac{q}{2}且最小系数严格大于q2-\frac{q}{2}即可,(经过上面计算可发现其实这两个是等价的,多了个负号而已),所以即需要保证:

(3d+12)p<q2(3d + \frac{1}{2}) · p < \frac{q}{2}

化简也即:

(6d+1)p<q(6d+1)·p < q

即需要在设置四个公开参数时保证以上关系,才能让a(x)a(x)正确地转换到RR上,保证解密的正确性。

接着按前文说的,现在a(x)a(x)RR上,可以转换到RpR_p上以消除pg(x)r(x)p · g(x) · r(x)

a(x)pg(x)r(x)+f(x)m(x)f(x)m(x) (in Rp)a(x) \equiv p · g(x) · r(x) + f(x) · m(x) \equiv f(x) · m(x) \ \text{(in $R_p$)}

为了恢复出m(x)m(x)还需要消去f(x)f(x),这直接乘上f1(x)f^{-1}(x)即可:

b(x)Fp(x)a(x)f1(x)f(x)m(x)m(x) (in Rp)b(x) \equiv F_p(x) · a(x) \equiv f^{-1}(x) · f(x) · m(x) \equiv m(x) \ \text{(in $R_p$)}

到这可以已经解密出m(x)m(x),但这个m(x)m(x)落在RpR_p,即系数取值在[0,p)[0, p),而实际的m(x)m(x)取值应该是在(p2,p2](-\frac{p}{2}, \frac{p}{2}],所以最后还要把b(x)b(x)做一次Center-lift才解出真正的m(x)m(x)

总结·

公开参数(N,p,q,d)(N, p, q, d),其中NN是素数、pp一般取33qq一般取22的幂,还需要保证gcd(N,q)=gcd(p,q)=1\rm{gcd}(N, q) = \rm{gcd}(p, q) = 1,和(6d+1)p<q(6d+1)·p < q

密钥生成:采样f(x)$T(d+1,d)f(x) \overset{\$}{\leftarrow} \mathcal{T}(d+1, d)g(x)$T(d,d)g(x) \overset{\$}{\leftarrow} \mathcal{T}(d, d),计算Fqf1(x) (in Rq)F_q \equiv f^{-1}(x) \ \text{(in }R_q \text{)}Fpf1(x) (in Rp)F_p \equiv f^{-1}(x) \ \text{(in }R_p \text{)},计算公钥h(x)Fq(x)g(x) (in Rq)h(x) \equiv F_q(x) g(x) \ \text{(in }R_q \text{)}

  • 私钥为:(f(x),g(x),Fq(x),Fp(x))(f(x), g(x), F_q(x), F_p(x))
  • 公钥为:h(x)h(x)

加密:明文m(x)m(x)的系数需落在(p2,p2](-\frac{p}{2}, \frac{p}{2}](虽然很奇怪下面表中说明文取自RpR_p),采样随机数r(x)$T(d,d)r(x) \overset{\$}{\leftarrow} \mathcal{T}(d, d),密文为e(x)ph(x)r(x)+m(x) (in Rq)e(x)\equiv ph(x)r(x) + m(x) \ \text{(in }R_q \text{)}

解密:计算a(x)f(x)e(x) (in Rq)a(x) \equiv f(x)e(x) \ \text{(in }R_q \text{)}并做Center-lift到a(x)Ra(x) \in R,计算m(x)Fp(x)a(x) (in Rq)m(x) \equiv F_p(x) a(x) \ \text{(in }R_q \text{)},最后把m(x)m(x)做Center-lift到RR(如果明文取自RpR_p则不需要这一步)。

参考代码·

参照上面理论部分写了个代码,其实我在很久以前就写过一个NTRU的代码,但因为穷拿去投题了,所以有保密协议就不能发,这次代码的加了很多骚操作,整体来说也更美观,所以感觉还是值得发一下。

首先是三个环RRRpR_pRqR_q,需要先做出Z[x]\mathbb{Z}[x](Z/pZ)[x](\mathbb{Z}/p\mathbb{Z})[x](Z/qZ)[x](\mathbb{Z}/q\mathbb{Z})[x],在Sage代码中写作:

1
2
3
R_  = PolynomialRing(ZZ,'x')
Rp_ = PolynomialRing(Zmod(p),'xp')
Rq_ = PolynomialRing(Zmod(q),'xq')

其中为了防止混淆,(Z/pZ)[x](\mathbb{Z}/p\mathbb{Z})[x](Z/qZ)[x](\mathbb{Z}/q\mathbb{Z})[x]的未知数在代码中我称作xpxq

然后使用上面三个东西做商环(模xN1x^N-1)就得到三个环RRRpR_pRqR_q

1
2
3
R  = R_.quotient(x^N - 1, 'y')
Rp = Rp_.quotient(xp^N - 1, 'yp')
Rq = Rq_.quotient(xq^N - 1, 'yq')

其中也是为了防混淆,我把RRRpR_pRqR_q上的未知数称作yypyq

然后采样三元多项式的函数T\mathcal{T},我选择先固定d1d_1个系数11d2d_2个系数1-1,然后再用shuffle将其打乱以获得随机的输出,然后要注意T\mathcal{T}的输出落在RR(毕竟有负数):

1
2
3
4
5
def T(self, d1, d2):
assert self.N >= d1+d2
t = [1]*d1 + [-1]*d2 + [0]*(self.N-d1-d2)
shuffle(t)
return self.R(t)

Center-lift函数基本照抄理论部分的公式就可以了,但由于作为函数我并不知道需要Lift的是RpR_p环还是RqR_q环,所以需要用“曲线救国”的方法在输入中提取模数(估计有更好的方法,懒):

1
2
3
def lift(self, fx):
mod = Integer(fx.base_ring()(-1)) + 1 # emmm
return self.R([Integer(x)-mod if x > mod//2 else x for x in list(fx)])

然后是密钥生成,首先f(x)f(x)g(x)g(x)直接用T\mathcal{T}采样即可。FpF_p的计算,由于在RpR_p上,而pp是素数,所以用Sage可以很方便地求逆:

1
Fp = self.Rp(list(fx)) ^ (-1)

上面需要先把fx转换到Rp上再做求逆,而由于前面设置的未知数名称不一致,不能直接用Rp(fx)做转换,曲线救国,先用list(fx)提取fx的系数,然后用系数转到Rp

FqF_q的计算则有点复杂,由于qq不为素数所以不能直接用Sage像上面那样求逆,目前网上查到的流行做法如这个教程invertmodpowerof2做递归求逆,但是我后来想到了一种更美观的“曲线救国”方法。

参考整数域上的求逆,比如要求a1(modn)a^{-1} \pmod n,除了用egcd求逆外,还有一种方法是(借助欧拉定理):

aφ(n)1aφ(n)a11a1a1(modn)a^{\varphi(n)-1} \equiv a^{\varphi(n)} a^{-1} \equiv 1 · a^{-1} \equiv a^{-1} \pmod n

其中的φ(n)\varphi(n)Zn\mathbb{Z}_n^*的阶。

所以只要求出FqF_q的阶即可利用类似方法求逆,问题是如何求FqF_q的阶。

经过计算加实验~~(原理未明,留坑)~~,发现f(x)f(x)FqF_q的阶(qq22的幂次)为:

φ(q)=qNq\varphi(q) = q^{N} - q

的因子。

PS:这里我取了最大值,即以上不是真正符合阶的定义,但是阶的某个倍数,所以说其阶是上面式子的一个因子,如果有方法可以分解上式则可能可以找到真正的阶,但在qq不大的情况下找到真正的阶其实并不能显著加速求逆的效率,所以可以选择直接使用φ(q)\varphi(q)

f(x)f(x)FpF_p的阶为(pp为素数):

φ(p)=pNp\varphi(p) = p^{N} - p

的因子。这个结果和之前求GL(n,Fp)GL(n, \mathbb{F}_p)阶的结果类似,毕竟pp为素数时,两个东西同构(如果我没记错的话,逃)。

PS:(2024.01.03更新)这里求阶的方法其实和RSA中求φ\varphi的方法类似,首先需要把所模的多项式分解为不可约的多项式,然后再分别求模这些不可约多项式的子环的阶,最后乘起来就是所求的真正的阶。但是这样算出来的阶会非常复杂,而pNpp^{N} - p又刚好阴差阳错地是这些阶的一个倍数,所以还是直接用pNpp^{N} - p为好(逃

另外,若ppqq不按[HPS14]的说明选择,根据中国剩余定理(CRT),假设p=ipikip = \sum_i p_i^{k_i},则f(x)f(x)FpF_p的阶为:

i(pikiNpiki)\prod_i (p_i^{k_iN} - p_i^{k_i})

综上,通过计算

Fpf(x)φ(p)1f1(x) (in Rp)Fqf(x)φ(q)1f1(x) (in Rq)F_p \equiv f(x)^{\varphi(p)-1} \equiv f^{-1}(x) \ \text{(in $R_p$)} \\ F_q \equiv f(x)^{\varphi(q)-1} \equiv f^{-1}(x) \ \text{(in $R_q$)}

即可算出FpF_pFqF_q

然后加密好像没啥可说的,只是为了统一,我把返回值转成字符串再输出,这样在输进解密函数时也不需要考虑类型问题。

解密函数也没啥好说,关注Center-lift的时机即可。

最后的参考代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# Sage

# Ref: https://www.osgeo.cn/sagemath/constructions/rings.html
class NTRU:
def __init__(self, N, p, q, d):
self.debug = False

assert q > (6*d+1)*p
assert is_prime(N)
assert gcd(N, q) == 1 and gcd(p, q) == 1
self.N = N
self.p = p
self.q = q
self.d = d

self.R_ = PolynomialRing(ZZ,'x')
self.Rp_ = PolynomialRing(Zmod(p),'xp')
self.Rq_ = PolynomialRing(Zmod(q),'xq')
x = self.R_.gen()
xp = self.Rp_.gen()
xq = self.Rq_.gen()
self.R = self.R_.quotient(x^N - 1, 'y')
self.Rp = self.Rp_.quotient(xp^N - 1, 'yp')
self.Rq = self.Rq_.quotient(xq^N - 1, 'yq')

# order check in keyGen
#self.RpOrder = self.p^(self.N - 1) - 1
#self.RqOrder = (self.q^self.N - self.q) // (self.q-1)
self.RpOrder = self.p^self.N - self.p
self.RqOrder = self.q^self.N - self.q

self.sk, self.pk = self.keyGen()


def test(self):
assert self.debug == True
pass


def T(self, d1, d2):
assert self.N >= d1+d2
t = [1]*d1 + [-1]*d2 + [0]*(self.N-d1-d2)
shuffle(t)
return self.R(t)

# center lift
def lift(self, fx):
mod = Integer(fx.base_ring()(-1)) + 1 # emmm
return self.R([Integer(x)-mod if x > mod//2 else x for x in list(fx)])


def keyGen(self):
fx = self.T(self.d+1, self.d)
gx = self.T(self.d, self.d)

Fp = self.Rp(list(fx)) ^ (-1) # list emmm
assert pow(self.Rp(list(fx)), self.RpOrder-1) == Fp # order checked
assert self.Rp(list(fx)) * Fp == 1

# Fq = self.Rq(fx) ^ (-1) # wasted
Fq = pow(self.Rq(list(fx)), self.RqOrder - 1) # invert
assert self.Rq(list(fx)) * Fq == 1 # order checked

hx = Fq * self.Rq(list(gx))

sk = (fx, gx, Fp, Fq, hx)
pk = hx
return sk, pk


def setKey(self, fx, gx):
assert type(fx) == type('x^2 + 1') # e.g.
assert type(gx) == type('x^2 - 1') # emmm

try:
fx = self.R(fx)
gx = self.R(gx)

Fp = self.Rp(list(fx)) ^ (-1)
Fq = pow(self.Rq(list(fx)), self.RqOrder - 1)
hx = Fq * self.Rq(list(gx))

self.sk = (fx, gx, Fp, Fq, hx)
self.pk = hx
return True
except:
return False


def getKey(self):
ssk = (
str(self.R_(list(self.sk[0]))), # fx
str(self.R_(list(self.sk[1]))) # gx
)
spk = str(self.Rq_(list(self.pk))) # hx
return ssk, spk


def encrypt(self, m):
assert type(m) == type('x^2 + 1') # e.g.
assert self.pk != None
hx = self.pk
mx = self.R(m)
mx = self.Rp(list(mx)) # change m to Rp, TODO: assert m in Rp
mx = self.Rq(list(mx)) # change m to Rq

rx = self.T(self.d, self.d)
rx = self.Rq(list(rx))


e = self.p * rx * hx + mx
#return e
return str(self.Rq_(list(e)))


def decrypt(self, e):
assert type(e) == type('xq^2 - 1') # e.g.
assert self.sk != None
fx, gx, Fp, Fq, hx = self.sk

e = self.Rq(e)
ax = self.Rq(list(fx)) * e
a = self.lift(ax) # center lift
bx = Fp * self.Rp(list(a))
b = self.lift(bx)

#return bx
return str(self.R_(list(b)))



if __name__ == '__main__':
mm = '-x^2 + x + 1'
ntru = NTRU(N=11, p=3, q=512, d=3)
#ntru.setKey('xp^2+1', 'xq^2-1')

print('keyGen check:')
sk, pk = ntru.getKey()
print("fx = '%s'" % sk[0])
print("gx = '%s'" % sk[1])
print("hx = '%s'" % pk)

print('\nencrypt/decrypt check:')
e = ntru.encrypt(mm)
print("e = '%s'" % e)
m = ntru.decrypt(e)
print(m)
assert m == mm
print(m)

print('\ncheck setKey:')
fx = 'x^10 + x^9 - x^8 - x^5 + x^4 + x - 1'
gx = '-x^10 + x^9 - x^6 - x^5 + x^4 + 1'
hx = '357*xq^10 + 156*xq^9 + 22*xq^8 + 467*xq^7 + 356*xq^6 + 23*xq^5 + 422*xq^4 + 490*xq^3 + 200*xq^2 + 201*xq + 378'
e = '286*xq^10 + 336*xq^9 + 355*xq^8 + 220*xq^7 + 330*xq^6 + 198*xq^5 + 182*xq^4 + 443*xq^3 + 454*xq^2 + 45*xq + 227'
ntru.setKey(fx, gx)
m = ntru.decrypt(e)
assert m == mm
print(m)

攻击方法·

下面说一下[HPS14]里提到的两种攻击方法。

小密钥空间·

[HPS14]的7.10.2中提到:

NNdd的取值不合理的话,会导致f(x)f(x)的取值空间过小,从而可以被攻击者枚举出密钥。

这个在参数设置上注意就好了,做题的话看见NNdd过小可以考虑直接枚举。

格规约攻击·

这里参考[HPS14]的7.11的内容,主要关注公钥的生成:

h(x)f(x)1g(x) in Rqh(x) \equiv f(x)^{-1} · g(x) \ \text{in } R_q

把模数去除,参考[HPS14]的写法即存在一个多项式u(x)u(x)使得:

f(x)h(x)qu(x)=g(x)f(x) · h(x) - q · u(x) = g(x)

写成矩阵形式也即(如果熟悉格规约的话,会把已知的放进矩阵):

(f(x),u(x))[1h(x)0q]=(f(x),g(x))(f(x), -u(x)) \begin{bmatrix} 1 & h(x) \\ 0 & q \\ \end{bmatrix} = (f(x), g(x))

这个只是一个写着方便的表述,实际上f(x)f(x)g(x)g(x)h(x)h(x)u(x)u(x)都是多项式,不能这样简单地占一个位置,应该

  • 多项式相乘:用前面说的多项式卷积乘法的矩阵写法;
  • 多项式数乘:用多项式的系数向量和对角线上为该常数的对角矩阵相乘;

这里我就直接抄书了:

用书中的简化写法是(为了混淆h(x)h(x)我用了HHII是单位矩阵):

MhNTRU=[IH0qI]M_h^{\rm{NTRU}} = \begin{bmatrix} \begin{array}{l|l} I & H \\ \hline 0 & qI \end{array} \end{bmatrix}

存在关系(记fff(x)f(x)的系数向量,记ggg(x)g(x)的系数向量):

(fu)MhNTRU=(fg)(f | -u) · M_h^{\rm{NTRU}} = (f | g)

直观上ffgg上元素的绝对值最大是11,所以(fg)(f | g)是“小”向量,用LLL对MhNTRUM_h^{\rm{NTRU}}规约即可出(fg)(f | g),这里先计算一下,我就直接借用[HPS14]的计算(懒):

但实际操作的情况是,LLL实际上解的是apprSVP,当矩阵规模变大时,求出(fg)(f | g)的能力会变小,所以只要设置足够大的NN即可防止这类攻击。

贴个攻击参考代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# Sage

def matrix_overview(BB):
for ii in range(BB.dimensions()[0]):
a = ('%02d ' % ii)
for jj in range(BB.dimensions()[1]):
if BB[ii, jj] == 0:
a += ' '
else:
a += 'X'
if BB.dimensions()[0] < 60:
a += ' '
print(a)

# cp NTRU.sage.py NTRU.py
from NTRU import NTRU

N = 11
q = 512
p = 3
d = 3
hx = '357*xq^10 + 156*xq^9 + 22*xq^8 + 467*xq^7 + 356*xq^6 + 23*xq^5 + 422*xq^4 + 490*xq^3 + 200*xq^2 + 201*xq + 378'
e = '286*xq^10 + 336*xq^9 + 355*xq^8 + 220*xq^7 + 330*xq^6 + 198*xq^5 + 182*xq^4 + 443*xq^3 + 454*xq^2 + 45*xq + 227'

Rq_ = PolynomialRing(Zmod(q),'xq')
h = vector(list(Rq_(hx)))
print(h)

M = matrix(ZZ, 2*N)
for i in range(N):
M[i, i] = 1
M[N+i, N+i] = q
for j in range(N):
M[i, N+j] = h[(j-i) % N]
matrix_overview(M)
#print(M)

L = M.LLL()
#L = M.BKZ(block_size=32)

def inT(f, d1, d2, N):
fl = list(f)
c1 = fl.count(1)
c_1 = fl.count(-1)
c0 = fl.count(0)
return c1==d1 and c_1==d2 and c1+c0+c_1 == N

R_ = PolynomialRing(ZZ,'x')
ntru = NTRU(N=N, p=p, q=q, d=d)
for i in range(N):
fg = L[i]
f = vector(ZZ, fg[:N])
g = vector(ZZ, fg[N:])

if inT(f, d+1, d, N):
pass
elif inT(-f, d+1, d, N):
f = -f
else:
continue

if inT(g, d, d, N):
pass
elif inT(-g, d, d, N):
g = -g # duoyu
else:
continue

print('Trying in i = %d' % i)
fx = str(R_(list(f)))
gx = str(R_(list(g)))
print('f = %s' % fx)
print('g = %s' % gx)
ntru.setKey(fx, gx)
m = ntru.decrypt(e)
print('m = %s' % m)

print('----------------\n')

有一个尴尬的情况是,MhNTRUM_h^{\rm{NTRU}}是一个N-Gap的格,也即LLL/BKZ规约后的矩阵中,前N行向量的规模是相近的,这会导致LLL/BKZ规约后的第一行不一定是(fg)(f | g),当然也会有(fg)(f | g)不是最短向量的情况。

曲线救国的方法是,把前N行向量都试一遍,先检测其是否落在T(d+1,d)\mathcal{T}(d+1, d),再尝试解密。

另外,实验测出有多个多个f(x)f(x)都可以实现解密的情况,合理怀疑一个公钥h(x)h(x)可以对应多个私钥f(x)f(x),也就是有可能上面解出来的前NN行都是对的,emmmm。