APP下载

一种新的面向硬件轻量级near-MDS矩阵的构造算法*

2019-09-10张怡帆任炯炯陈少真

密码学报 2019年4期
关键词:分支线性乘法

张怡帆,任炯炯,陈少真

1.解放军信息工程大学,郑州450001

2.数学工程与先进计算国家重点实验室,郑州450001

1 引言

线性扩散层是对称密码的一个重要组件,它为对称密码算法提供了内部独立性.对于对称密码算法来说,其线性扩散层是设计时需要重点考虑的关键部件.我们用分支数来衡量一个扩散层的安全性,分支数较大的扩散层能更好地抵抗差分和线性攻击.在资源受限的环境下,不仅需要考虑安全性,而且要关注轻量级密码算法扩散层的硬件效率.随着轻量级密码的快速发展[1–4],如何构造分支数较大,硬件效率高的轻量级扩散层引起人们的重视.

分支数达到最大的矩阵称作MDS矩阵,能实现最好的扩散效果,因此许多MDS矩阵如递归MDS矩阵和对合MDS矩阵的构造方法被提出[5–16].许多算法诸如AES 等的线性扩散层采用MDS矩阵以期达到最好的扩散性,但MDS矩阵中的每一个元素均为非零,因此在硬件应用时耗能较大,效率较低.在资源受限的今天,算法实现起来较为困难.随着物联网技术等普适计算的发展,对硬件提出了新的安全需求,轻量级密码的研究领域越来越受重视.轻量级密码算法对安全性的要求相较来说有所降低,而如何在保证安全性的前提下提高效率受到了广泛关注.

2008年,Choy[17]和Khoo[9]将差分分支数达到次优的矩阵定义为almost-MDS矩阵,由于两篇文章中的almost-MDS矩阵对应不同文献给出的线性码定义,因此将差分与线性分支数均达到次优的矩阵称为near-MDS矩阵.与MDS矩阵相比,near-MDS矩阵的分支数仅次于MDS矩阵,安全性也仅次于MDS矩阵,但其具有实现代价小的优势.在对安全性的要求降低,实现代价要求提高的资源受限环境下,用near-MDS矩阵作为扩散层是一个很好的选择.实际上,一些用near-MDS矩阵构造的扩散层在硬件应用上的表现优于MDS矩阵或者迭代MDS矩阵.近来,near-MDS矩阵已经被许多轻量级分组密码采用,包括PRINCE[18]、FIDES[19]、PRIDE[20]、Midori[21]和MANTIS[22]等.Near-MDS矩阵被广泛应用在低消耗低延迟的轻量级分组密码的设计中,因此研究此类矩阵具有一定的密码学价值和现实意义.

有限域F2m中的每个元素都可以表示成F2上的m×m矩阵,所以F2m上的MDS矩阵实际上对应着F2上的一个分块矩阵,每一块都是m×m大小,需要注意的是,并不是F2上的每一个矩阵都可以代表有限域F2m内的一个元素,这取决于它的最小多项式是否为不可约的.Li[23]利用有限域F2m内的元素作为扩散矩阵的元素,在F24上构造出10 个XORs 最少的4×4 near-MDS矩阵.但是有限域F2m内的元素乘法问题[5]只是上线性变换的一种特殊形式,在上还存在着许多不能用有限域F2m中的元素乘法表示的线性变换[10],上的线性变换空间可以表示为F2上的m×m矩阵.因此,若只考虑在有限域F2m中构造MDS矩阵有可能会导致丢失一些较好的构造结果,在上的线性变换空间中,有可能构造出新的符合要求的轻量级near-MDS矩阵.在之前的构造方法中,用来构造矩阵的元素满足成对交换,或者假定满足成对交换[24,25],但该假定条件也可能会导致丢失较好的构造结果,所以本文中,没有假定用来构造near-MDS矩阵的上的线性变换是成对交换的.

在对称密码算法中,最常用的是4 比特S盒和8 比特S盒,最常用扩散矩阵的维度是4.当矩阵维度较大时,计算复杂度会变得非常复杂,因此将重点放在利用和上的线性变换来构造4×4 near-MDS矩阵上.本文利用直接将F2上的m×m矩阵作为扩散矩阵元素的方法,利用循环矩阵、对合矩阵的性质给出了循环对合矩阵元素间满足的性质引理,利用引理给出算法搜索的约束条件,并借助Matlab 构造出4×4 循环对合near-MDS矩阵,在m=4时可以找到48 个满足条件XORs 最少的4×4 循环对合near-MDS矩阵,m=8时,依据文献[9,14,26]中提到的子域构造方法,可直接得到异或操作数最少的8×8 循环对合near-MDS矩阵,较之前结果找到了更多的XORs 最少的循环对合near-MDS矩阵,且计算复杂度较小,在Windows 10 系统、i5-6200U CPU 处理器、4 G 内存的机器条件下仅需要大概6 分钟.

本文结构安排如下:第2 节给出了本文所涉及的符号表示及相关基础知识; 第3 节给出对near-MDS矩阵内元素的一个基本限制条件以及对合循环near-MDS矩阵的构造方法,并给出了自动搜索算法以及结果; 第4 节给出了一个简短的总结.附录中记录了本文所有的搜索结果.

2 基础知识

2.1 符号表示

XORs:异或操作次数;

Bd(L):矩阵L的差分分支数;

Bl(L):矩阵L的线性分支数;

#A:A·x所需的异或操作次数;

ω(A[i]):矩阵A中第i行非零元素的个数.

2.2 相关基础知识

本节介绍一些相关的基础概念和定义.

2.2.1 分支数

对x,y∈,若映射A:→满足A(x+y)=A(x)+A(y),则称为线性映射.在F2上固定的一组基,则上的一个线性映射可以表示为F2上的一个m×m矩阵,我们用A来表示,满足A(x)=A·x,其中x=(x1,··· ,xm)∈,在本文中代表一个列向量.一个线性映射是上的一个置换当且仅当它的矩阵表示是非奇异的.GL(m,S)代表所有m×m的非奇异矩阵,矩阵中元素属于S,S是一个元素集.

每个线性扩散层都是一个线性映射,可以用如下矩阵表示:

其中Li,j是F2上的m×m矩阵,对X=(x1,··· ,xn)∈L(X)=其中对于X=(x1,···,xn)∈中非零元素的个数定义为重量,用

定义1矩阵L是上的n×n矩阵,L的差分分支数定义为

L的线性分支数定义为

分支数可以用线性码的最小距离来表示.

定理1[27]假设矩阵L是上的n×n矩阵,C是上的一个[2n,n]线性码,它的生成矩阵是(In|LT),其中In是一个n维单位矩阵.则矩阵L的差分分支数等价于线性码C的最小距离,即Bd(L)=d(C).此外,Bl(L)=d(C⊥),其中C⊥是C的共轭码.

令C是一个[N,k]线性码,当C的最小距离达到上界d(C)=N−k+1[28]时,我们称其为一个MDS 线性码.若生成矩阵是(Ik|LT)的线性码CL是一个MDS 线性码,则称矩阵L是MDS矩阵.MDS矩阵的差分分支数和线性分支数同时达到上界[29],即Bl(L)=d(CL)=N−k+1.

MDS矩阵分支数达到最大,能实现最好的扩散效果,但是由于MDS矩阵的每一个元素均为非零,在硬件应用时耗能较大,效率较低.为折衷考虑安全性与效率的问题,本文主要研究非MDS矩阵中分支数达到最大的矩阵,下面我们给出near-MDS矩阵的定义.

定义2当Bd(L)=Bl(L)=d(CL)=n时,上的n×n矩阵L是near-MDS矩阵.

在文献[29]中,一个[n,k]线性码C满足条件d(C)=n−k,d(C⊥)=k.那么根据定理1,若一个n×n矩阵L满足Bd(L)=Bl(L)=n,那么矩阵(In|LT)是一个[2n,n,n]near-MDS 线性码的生成矩阵.由此我们给出定理2.

定理2[30]L是一个n维的非MDS矩阵,n是一个正整数且n≥2.则L是一个near-MDS矩阵的充要条件是:对任意的1≤g≤n−1,矩阵L的每一个g×(g+1)子矩阵和(g+1)×g子矩阵都至少有一个非奇异的g×g子矩阵.

由定理2 可知,当n较大时,为验证near-MDS矩阵的性质,计算量会相当复杂.因此本文重点研究4×4 矩阵,这是在密码算法中被广泛使用的矩阵大小.本文利用循环矩阵来构造轻量级near-MDS矩阵,循环矩阵可以用第一行的元素来定义,因此可以重复利用元件,提高硬件应用效率.

定理3[11]对任意置换矩阵P1和P2,矩阵L和P1LP2有相同的差分和线性分支数.

由定理3 可得,当循环矩阵Circ(A,B,C,D)的第一行元素发生位置改变时,矩阵的分支数保持不变.

2.2.2 XOR 数

a,b∈F2,a+b被称为比特XOR 运算.对于A∈GF(m,F2),#A表示A·x所需的XOR 操作次数,其中表示矩阵A中第i行非零元素的个数.例如为方便描述,本文用矩阵每一行非零元素的位置给出一个矩阵的简单表达式,例如[2,3,4,[1,4]]是矩阵A的简单表达式.

对之前给出的矩阵L=令X=(x1,··· ,xn)∈我们给出所需要的异或操作数,其中Li,j(xj)=以矩阵L的第一行为例,所需的异或操作数为#L1,1+#L1,2+···+#L1,n+(l−1)·m,l为第一行中非零Li,j的个数.

2.2.3 循环矩阵

特殊矩阵具有特殊的结构特点,其特点在应用过程中可能有利于降低硬件耗能,因此考虑借助特殊矩阵形式来构造扩散矩阵.本小节给出本文所用特殊矩阵的定义.

定义3一个矩阵被称为循环矩阵当且仅当它的每一行都是前一行向右(向左)循环一个元素的位置.

例如对一个4×4 的循环矩阵我们表示为:Circ(A,B,C,D)=循环矩阵每一行的元素均与第一行相同,在硬件应用时可以重复利用乘法元件,进而大大提高硬件应用效率,因此被很多轻量级密码算法所采用.

3 轻量级循环near-MDS矩阵

本节研究轻量级循环对合near-MDS矩阵的构造方法,证明矩阵元素满足循环对合性质所必须遵循的性质引理,根据引理,给出了搜索算法的约束条件条件,借助Matlab 寻找满足条件的异或数最少的矩阵.本文搜索到了48 个满足条件的循环对合near-MDS矩阵,较之前搜索结果更多,并且减少了搜索的时间复杂度,在Windows 10 系统、i5-6200U CPU 处理器、4 G 内存的机器条件下需要大概6 分钟.

3.1 Near-MDS矩阵的元素限制条件

near-MDS矩阵的分支数达到次优,其矩阵元素有如下特点:

引理1假设循环矩阵Circ(A,B,C,D)是一个near-MDS矩阵,则矩阵元素A,B,C,D中至多有一个为零矩阵O.

证明:若A,B,C,D中有两个元素为O,假设为A,B.根据定理2,当g=1时矩阵Circ(A,B,C,D)有一个1×2 的子矩阵不满足至少有一个1×1 的子矩阵是非奇异的,与定理2 给出的near-MDS 性质矛盾.所以矩阵元素A,B,C,D中至多有一个为零矩阵O.

本文的主要目标是减少搜索范围,寻找 XORs 最少的 near-MDS矩阵.根据引理1,假设Circ(A,B,C,D)中一直存在一个零矩阵O,根据定理3,循环矩阵第一行元素位置发生变化时,矩阵的分支数保持不变,因此本文重点研究Circ(O,A,B,C)这种形式的矩阵.

3.2 构造循环对合矩阵

循环矩阵元素排列有很明显的特点,若要一个循环矩阵Circ(O,A,B,C)同时满足对合性,其元素相互之间也存在一定的约束条件,对此本文有如下结果:

引理2L=Circ(O,A,B,C)是一个循环矩阵,其中A,B,C∈GL(m,F2).矩阵L满足对合性当且仅当下列等式同时满足:AB=BA,BC=CB,A2=C2,AC+CA+B2=I.

证明:通过矩阵乘法,可以得到

若矩阵L具有对合性,则满足条件L2=Circ(I,0,0,0),因此L是一个对合矩阵当且仅当下列条件同时满足:AB=BA,BC=CB,A2=C2,AC+CA+B2=I.

由引理2 给出的矩阵Circ(O,A,B,C)满足对合性其元素之间具有的约束条件,下面给出循环对合矩阵的一般构造方法.首先我们给出乘法阶的定义:对于矩阵A∈GL(m,F2),满足等式Ad=I的最小正整数d称为A的乘法阶.

引理3假设A,C∈GL(m,F2)且满足条件A2=C2=I,I+A+C的乘法阶为4k−2,其中k≥1.令B=(I+A+C)2k,则矩阵Circ(O,A,B,C)是一个对合矩阵.

证明:令B=(I+A+C)2k,因为A和C满足条件A2=C2=I,根据引理2,我们只需证明下面等式即可:

首先容易证明

然后可得B=(I+A+C)2k=(I+AC+CA)k.因此

同理我们可以证得BC=CB.因为(I+A+C)4k−2=I,因此我们可以得到

根据引理3,我们得到了矩阵Circ(O,A,(I+A+C)2k,C)是对合矩阵.

注意:当A=C时,(I+A+C)的乘法阶为1,不满足4k−2 的条件.这种情况下L=Circ(O,A,B,C)=Circ(O,A,B,A)仍为一个循环矩阵,若要满足对合性应同时满足条件B2=I,AB=BA,根据条件仍可寻找循环对合矩阵.

由引理2、引理3,我们给出要使矩阵满足循环对合性,矩阵元素所需满足的条件,利用引理给出自动搜索算法的约束条件,进而可以很大程度减小搜索范围.再利用定理2 给出的near-MDS矩阵的特性,借助Matlab,进一步搜索循环对合near-MDS矩阵.下面给出算法介绍以及搜索结果.

3.3 自动化搜索算法及结果

在上一节中我们推导出满足循环对合条件的矩阵元素之间所需满足的条件

(1)A=C时,A2=C2=I,(I+A+C)4k−2=I,k≥1,B=(I+A+C)2k.

(2)A=C时,B2=I,AB=BA.

将矩阵元素所满足的式子转化为算法中的约束条件,缩小搜索范围,在得到循环对合矩阵之后进一步检验其是否具有near-MDS矩阵的性质,最终获得循环对合near-MDS矩阵.构造算法步骤如下:

首先获得一个集合S,其中包含所有满足对合性的矩阵.然后分以下两种情况讨论:

(1)A=C时,我们选取矩阵对(A,C)∈S×S,计算I+A+C的乘法阶d,若dmod 4=2,则是循环对合矩阵.

(2)A=C时,选取B∈S,若满足AB=BA,则Circ(O,A,B,A)是循环对合矩阵.

本文借助Matlab 计算机程序,通过添加对矩阵元素的限制条件缩小搜索范围,对满足条件的循环对合near-MDS矩阵进行算法搜索.为寻找异或数最少的矩阵,矩阵元素集S的选择在搜索算法中起到了关键的作用,一般来说,元素的XORs 越少,矩阵的XORs 也会相应较少,因此在选择矩阵元素集时,首先从XORs 最少的矩阵集S开始搜索.因为XORs 最少为0,所以本文首先寻找满足对合性且XORs 为0即非零元素个数为4 的矩阵,将其作为矩阵集S,在S中搜索是否存在满足上述两种情况之一(A,B,C).若不存在,则扩大矩阵集S的范围,寻找满足对合性且XORs 为0 或者XORs 为1 的矩阵,将其作为矩阵集S进行搜索.以此类推,直到找到满足条件的(A,B,C).

其次,在得到循环对合矩阵后我们要验证其是否为near-MDS矩阵,根据定理2 给出的near-MDS 性质,对任意的1≤g≤n−1,矩阵L的每一个g×(g+1)子矩阵和(g+1)×g子矩阵都至少有一个非奇异的g×g子矩阵.若一个矩阵为非奇异矩阵,则其行列式一定不为0.为验证near-MDS 性质,对于一个g×(g+1)子矩阵或者(g+1)×g子矩阵,我们需要计算其所有的g×g子矩阵行列式的绝对值,为保证至少有一个非奇异的g×g子矩阵,应满足所有的g×g子矩阵的行列式绝对值之和大于0.由此,为验证所构造循环对合矩阵是否为near-MDS矩阵,本节给出如下引理及证明.

引理4令d1,d2,··· ,dk≥0 代表k个正数值,其中k是一个正整数,若d1,d2,··· ,dk中至少有一个非零值当且仅当d1+d2+···+dk>0.

证明:因为d1,d2,··· ,dk≥0,要满足d1+d2 +···+dk=0,当且仅当d1=d2=···=dk=0,与所给条件矛盾.因此应满足d1+d2+···+dk>0.

综合考虑之前给出的元素约束条件以及near-MDS矩阵性质验证方法,利用Matlab 计算机程序,给出的具体搜索算法.见算法1.

算法1 构造循环对合矩阵并检验near-MDS矩阵的性质Input:矩阵元素限制条件; 4×4 循环矩阵Circ(O,A,B,C);Output:循环对合near-MDS矩阵Circ(O,A,B,C)1 S ←∅;2 if N2=I & & #N=4 then 24 for g ∈[1,n−1]do 25 for 矩阵T 所有g×(g+1)和(g+1)×g 子矩阵M do 3 N ∈S;4 end 5 for A=C do 26 p ←0;27 for 矩阵M 的所有g×g 子矩阵N do 6(A,C)∈S×S;7 计算I +A+C 的乘法阶d ;8 if d mod 4=2 then 28 计算行列式det(N)的绝对值q;29 p ←p+q;30 if p=0 then 9 B=(I +A+C)d2+1;10 end 11 end 12 for A=C do 31 return false;32 end 33 else 13 B ∈S;14 if AB=BA then 15 T=Circ(O,A,B,C);16 end 17 end 18 if 不存在(A,B,C)满足条件then 34 T=Circ(O,A,B,C)是一个循环对合35 near-MDS矩阵36 end 37 end 38 end 39 end 40 return T 19 重新选择矩阵集S;20 if N2=I & & #N=4+1=5 then 21 N ∈S;22 end 23 end

搜索结果如下:

m=4时,本文选取S∈GL(4,F2).首先搜索XORs 为0 的矩阵集S,在S中搜索到了满足Circ(O,A,B,C)是near-MDS矩阵的(A,B,C),并且#A+#B+#C=0,达到了矩阵第一行元素异或操作数之和的最小值.利用上述算法我们找到了48 组满足条件的(A,B,C),其中有36 组满足A=C,12 组满足A=C.

实例:

(1)m=4 且A=C时,A=[4,2,3,1],C=[4,2,3,1],B=[1,3,2,4].

(2)m=4 且A=C时,A=[3,4,1,2],C=[4,3,2,1],B=[1,2,3,4].

m=8时,因为已经构造出m=4时满足#A+#B+#C=0 的循环对合near-MDS矩阵,根据文献[9,14,26]中提到的子域构造方法,我们很容易在GL(8,F2)中构造出满足#A+#B+#C=0 的循环对合near-MDS矩阵.举例介绍子域构造方法,令X,Y,Z∈GL(4,F2),#X=#Y=#Z=0 并且Circ(O,X,Y,Z)是一个循环对合near-MDS矩阵,则对A,B,C∈GL(8,F2),Circ(O,A,B,C)也是一个循环对合near-MDS矩阵,其中

实例:

(1)m=8 且A=C时,A=[4,2,3,1,8,6,7,5],C=[4,2,3,1,8,6,7,5],B=[1,3,2,4,5,7,6,8].

(2)m=8 且A=C时,A=[3,4,1,2,7,8,5,6],C=[4,3,2,1,8,7,6,5],B=[1,2,3,4,5,6,7,8].

表1 给出了用本文方法构造的满足条件的矩阵数与已知的用有限域内元素乘法问题构造的矩阵数的比较,可以看出本文构造出了较之前结果更多的满足条件的矩阵,这对算法的设计具有一定的借鉴意义.所有的构造算法结果见附录.

表1 本文结果与文献[22]结果对比Table 1 Compare results of proposed method and Ref.[22]

4 总结

在资源受限的环境下,轻量级密码算法扩散层的设计更加注重硬件效率.随着轻量级密码的快速发展,构造分支数较大,硬件效率高的轻量级扩散层是一个亟需解决的问题.

本文主要探索了用F2上的m×m矩阵来构造4×4 轻量级near-MDS矩阵的方法,理论上证明了如何判断循环对合的性质引理,在此基础上,将引理条件转化为自动搜索算法的约束条件,给出了Matlab程序编译算法,实现了对满足条件矩阵的算法搜索.利用循环矩阵和对合矩阵的特性推导出矩阵元素之间的关系,在Matlab 算法搜索过程中添加矩阵元素的约束条件,缩小搜索范围,降低复杂度.同时选择XORs 较少的矩阵元素,以获得硬件效率高的扩散矩阵.本文对F2上的4×4 循环对合near-MDS矩阵进行了搜索,找到了48 个XORs 最少的矩阵,较之前利用有限域内元素乘法构造的满足条件的10 个矩阵,不仅给出满足条件的矩阵数量更多,而且降低了搜索的时间复杂度.如何利用其它类型特殊矩阵的性质获得矩阵元素约束条件以及构造更大维度的轻量级near-MDS矩阵还需要进一步的研究与实验.

附录

(1)A=C,m=4时,满足B2=I、AB=BA,Circ(O,A,B,C)是上的循环对合near-MDS矩阵.我们找到了36 组满足条件的[A,B].

1.A=[4,3,2,1]B=[4,3,2,1]

2.A=[4,3,2,1]B=[3,4,1,2]

3.A=[4,3,2,1]B=[2,1,4,3]

4.A=[4,3,2,1]B=[1,2,3,4]

5.A=[4,2,3,1]B=[4,2,3,1]

6.A=[4,2,3,1]B=[1,3,2,4]

7.A=[3,2,4,1]B=[1,2,3,4]

8.A=[2,4,3,1]B=[1,2,3,4]

9.A=[3,4,1,2]B=[4,3,2,1]

10.A=[3,4,1,2]B=[3,4,1,2]

11.A=[3,4,1,2]B=[2,1,4,3]

12.A=[3,4,1,2]B=[1,2,3,4]

13.A=[4,2,1,3]B=[1,2,3,4]

14.A=[3,2,1,4]B=[3,2,1,4]

15.A=[3,2,1,4]B=[1,4,3,2]

16.A=[2,3,1,4]B=[1,2,3,4]

17.A=[4,1,3,2]B=[1,2,3,4]

18.A=[3,1,2,4]B=[1,2,3,4]

19.A=[2,1,4,3]B=[4,3,2,1]

20.A=[2,1,4,3]B=[3,4,1,2]

21.A=[2,1,4,3]B=[2,1,4,3]

22.A=[2,1,4,3]B=[1,2,3,4]

23.A=[2,1,3,4]B=[2,1,3,4]

24.A=[2,1,3,4]B=[1,2,4,3]

25.A=[1,4,3,2]B=[3,2,1,4]

26.A=[1,4,3,2]B=[1,4,3,2]

27.A=[1,3,4,2]B=[1,2,3,4]

28.A=[1,4,2,3]B=[1,2,3,4]

29.A=[1,3,2,4]B=[4,2,3,1]

30.A=[1,3,2,4]B=[1,3,2,4]

31.A=[1,2,4,3]B=[2,1,3,4]

32.A=[1,2,4,3]B=[1,2,4,3]

33.A=[1,2,3,4]B=[4,3,2,1]

34.A=[1,2,3,4]B=[3,4,1,2]

35.A=[1,2,3,4]B=[2,1,4,3]

36.A=[1,2,3,4]B=[1,2,3,4]

(2)A=C,m=4时,(I+A+C)2=I,此时B=(I+A+C)2=I,Circ(O,A,B,C)是上的循环对合near-MDS矩阵.我们找到了12 组满足条件的[A,C].

1.A=[4,3,2,1]C=[3,4,1,2]

2.A=[4,3,2,1]C=[2,1,4,3]

3.A=[4,3,2,1]C=[1,2,3,4]

4.A=[3,4,1,2]C=[4,3,2,1]

5.A=[3,4,1,2]C=[2,1,4,3]

6.A=[3,4,1,2]C=[1,2,3,4]

7.A=[2,1,4,3]C=[4,3,2,1]

8.A=[2,1,4,3]C=[3,4,1,2]

9.A=[2,1,4,3]C=[1,2,3,4]

10.A=[1,2,3,4]C=[4,3,2,1]

11.A=[1,2,3,4]C=[3,4,1,2]

12.A=[1,2,3,4]C=[2,1,4,3]

猜你喜欢

分支线性乘法
算乘法
一类离散时间反馈控制系统Hopf分支研究
软件多分支开发代码漏合问题及解决途径①
二阶整线性递归数列的性质及应用
《整式的乘法与因式分解》巩固练习
把加法变成乘法
巧分支与枝
非齐次线性微分方程的常数变易法
硕果累累