CN108228138A - 一种sidh中特殊域快速模乘的方法 - Google Patents

Abstract

本发明提出一种SIDH中特殊域快速模乘的方法，该方法对现有的SIDH中特殊域快速模乘计算方法进行了改进，通过构建关系式：将SIDH中单个模乘A·B(mod p)的计算转换为A′·B′(mod p)的计算，在计算过程中大大节省了操作数的数目，从而提高计算速度。

Description

一种SIDH中特殊域快速模乘的方法

技术领域

本发明涉及SIDH中特殊域模乘计算领域，尤其是一种SIDH中特殊域快速模乘的方法。

背景技术

在广泛的后量子技术之中，最新的超奇异同源点Diffie-Hellman(SIDH)密钥交换协议具很多优良特性以及广泛的应用前景。与其他广泛被引用的后量子密钥交换和加密体系相比，如格加密代码加密、哈希加密以及多变量加密，SIDH密钥交换体系所需要的 密钥长度要小得多。2011年，Jao和Deo提出了一种基于SIDH中特殊结构域的快速模 乘算法，该算法中的一个计算瓶颈就是p＝f·2^a3^b-1的模乘运算，在该算法中很难兼 顾资源消耗量和运算吞吐量之间的平衡，即模乘运算所占用的资源较多，但运算吞吐量 并不大。

发明内容

发明目的：为解决上述技术问题，本发明提出一种SIDH中特殊域快速模乘的方法。

技术方案：为解决上述技术问题，在增加较少资源占用率的前提下获得更快的处理 速度，本发明提出一种SIDH中特殊域快速模乘的方法，该方法包括SIDH中单个模乘的计算步骤：

(1)定义SIDH中单个模乘为A·B(modp)，其中，p为素数域，A和B为素数域 p中的任意两个元素；p＝2·2^a3^b-1，a，b均为偶数；

(2)将A·B(modp)表示为：

式中，A、A′、B、B＇分别用以R为基的补码形式表示为：

A＝a₁R²+a₂R+a₃

B＝b₁R²+b₂R+b₃

A′＝a＇₂R+a′₃

B′＝b′₂R+b′₃

其中，R＝2^a/23^b/2，a₁、a₂、a₃、b₁、b₂、b₃、a＇₂、a′₃、b′₂、b′₃均为系数，符号表示异或；

(3)通过计算A′·B′(modp)得到A·B(modp)，包括依次执行步骤(3-1)至(3-8)：

(3-1)将A′·B′(modp)的值以基R形式表示为：

A′·B′(modp)＝c′₁R²+c′₂R+c′₃

式中，c＇₁、c′₂、c′₃为系数变量；

(3-2)初始化c＇₁＝0，c′₂＝0，c′₃＝0；

(3-3)计算c＇₁＝c′₂c′₃[0]，[0]表示取最低位；

(3-4)计算

(3-5)对c′₃进行Barrett约减，约减结果记为r₁，即c′₃＝r₁；计算c′₂＝c′₂+r₁；

(3-6)对c′₂进行Barrett约减，约减结果记为r₂，即c′₂＝r₂；计算c＇₁＝c＇₁+r₂；

(3-7)计算c＇₁＝c＇₁[0]；

(3-8)根据计算得到的c＇₁、c′₂、c′₃得到A′·B′(modp)的计算结果，再根据 A′·B′(modp)得到A·B(modp)的计算结果。

有益效果：与现有技术相比，本发明具有以下优势：

本发明可以节省操作数的数目，从而提高计算速度。在硬件实现时可以在增加较少 硬件资源的情况下，将SIDH的吞吐量提高到了原来的6倍以上。

附图说明

图1为本发明的原理图；

图2为Barrett约减算法原理图；

图3为Barrett除法算法原理图；

图4为现有的SIDH中特殊结构域的模乘计算原理图；

图5为本发明的硬件架构图；

图6为本发明中乘法操作的硬件结构图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

一、现有技术说明

图4所示为现有的SIDH中特殊结构域的模乘计算原理图，下面首先介绍现有的SIDH

中特殊结构域的模乘计算原理。

假定p＝2·2^a3^b-1，其中a和b为偶数，通过基R表示法，可以将素数域p中的任 一元素(即所有小于p的非负整数)A表示成如下形式，其中R＝2^a/23^b/2，

A＝a₁R²+a₂R+a₃，a₁∈{0，1}，a₂，a₃∈[0，R) (1)

式中，a₁、a₂、a₃均为系数。

将p中的任意两个元素A和B均通过基R表示，那么A·B(modp)的结果为：

C＝a₁b₁·R⁴+(a₁b₂+a₂b₁)·R³+(a₁b₃+a₂b₂+a₃b₁)·R²+(a₂b₃+a₃b₂)·R+a₃b₃ (2)

由于p＝2·2^a3^b-1，可以得到公式(2)中2^a3^b亦即R²＝2^-1(mod p)，因此a₁b₁·R⁴或是等于0或是等于2^-2(mod p)，因为a₁，b₁∈{0，1}所以a₁b₁∈{0，1}。当选定某个固定不变 的p之后，对应的2^-2(mod p)的值也是固定不变的，在算法开始之前就预计算出 2^-2(mod p)的值，从而减少算法中的操作。

同样，对于其他项，也可以用类似的方法代替。例如，对于(a₁b₃+a₂b₂+a₃b₁)·R²项，如果该项是偶数，可以得到：

(a₁b₃+a₂b₂+a₃b₁)·R²＝(a₁b₃+a₂b₂+a₃b₁)/2(mod p)

反之，若为奇数，则可以写成

(a₁b₃+a₂b₂+a₃b₁)·R²＝(a₁b₃+a₂b₂+a₃b₁-1)/2(mod p)+((a₁b₃+a₂b₂+a₃b₁)mod 2)·2^a3^b

将两种情况结合起来考虑便可以得到：

其中，符号表示向下取整。

类似的，对于(a₁b₂+a₂b₁)·R³项，可以写成

因此，公式(2)可以重新写成如下形式：

但是在式(3)之中，最后得到的系数c₂和c₃是在[0，R²)范围之内，而不是在基R示法中 要求的[0，R)范围之内，因此，需要对c₂和c₃作进一步化简。

对c2和c3的化简采用了Barrett约减算法，Barrett约减算法的原理如图2所示：

根据欧拉除法引理可得，对于任意两个正整数a和b，必定存在q和r使得等式 a＝q·b+r成立，其中r∈[0，b-1]，亦即a＝r(mod b)。

毫无疑问，得到这样的q和r必然需进行一次a/b的除法，但是在计算机中，相比 于乘法，除法是代价很大的操作。然而对于常除数而言，Barrett约减却是一种很巧妙的 操作。它可以将一次1/b的除法转化成若干乘法和移位操作。在Barrett约减中，1/b可 以表示成：

通常来说，x的取值为其中参数k的取值取决于a。从公式(4)中可以 看出，对1/b作近似化后产生的误差值e＝1/b-x/2^k。因此，对于商q产生的总误差值 为ae。由于q∈Z⁺，为了使最终得到的结果正确，需要使误差ae小于1，因此需要满足 条件k＝log₂(a)。

在对c₂和c₃进行化简时，需要对c₂和c₃做关于R的除法操作，得到相对应的商q 和余数r。由于R＝2^a/22^b/2，可以将该除法看成一个先除以2^a/2之后再除以3^b/2的操作。 而在计算机中，关于2^a/2的除法可以用简单的右移操作代替。因此对c₂和c₃的快速化简， 可以按照以下步骤进行：

1、提取出c_i的低a/2位并将它存在变量r₁中；

2、对c_i右移a/2位得到c′_i；

3、对c′_i作除以3^b/2的除法得到商q以及余数r₂。

4、得到c_i＝q·2^a/23^b/2+(r₂·2^a/2+r₂)＝q·2^a/23^b/2+r

然而，步骤3中的关于3^b/2的除法操作并不像关于2^a/2的除法操作一样能简单实现，由于被除数3^b/2一直是固定不变的，因此考虑用Barrett约减算法来完成关于3^b/2的除法 操作部分。综上所述，对于c₂和c₃的完整化简过程，称之为Barrett除法算法，Barrett 除法算法的过程如图3所示。

至此，我们可以得到完成SIDH中一次模乘A·B(modp)的计算过程。

二、本发明技术方案的说明

根据模运算的基本性质，可以得到以下等式：

(p-A)·(p-B)(mod p)＝A·B(mod p) (5)

首先，同样设任一元素A∈F_p，F_p为素数域；把A以基R形式表示为：

A＝a₁R²+a₂R+a₃，a₁∈{0，1}，a₂，a₃∈[0，R)

可以得到元素A′∈F_p：

当a₁＝1时，有

a′_i＝R-a_i，i∈{2，3} (7)

则A′的基R形式表达式为：

A′＝a‘₂R+a’₃，a₂，a₃∈[0，R) (8)

同理，得到B′。

根据公式(5)可以得到算式A·B(modp)和A′·B′(modp)满足以下关系：

其中符号表示异或。

通过计算A′·B′(mod p)即可得到A·B(mod p)，计算A′·B′(mod p)的步骤为：

(1)将A′·B′(mod p)的值以基R形式表示为：

A′·B′(modp)＝c′₁R²+c′₂R+c′₃

式中，c＇₁、c′₂、c′₃为系数变量；

(2)初始化c＇₁＝0，c′₂＝0，c′₃＝0；

(3)计算c＇₁＝c′₂c′₃[0]，[0]表示取最低位；

(4)计算

(5)对c′₃进行Barrett约减，约减结果记为r₁，即c′₃＝r₁；计算c′₂＝c′₂+r₁；

(6)对c′₂进行Barrett约减，约减结果记为r₂，即c′₂＝r₂；计算c＇₁＝c′₁+r₂；

(7)计算c′₁＝c＇₁[0]；

(8)根据计算得到的c′₁、c′₂、c′₃得到A′·B′(mod p)的计算结果，再根据A′·B′(mod p) 得到A·B(modp)的计算结果。

与现有的SIDH中一次模乘计算方法相比，本发明和现有技术都需要计算4个乘法：原算法的a₂×b₂，a₂×b₃，a₃×b₂以及a₃×b₃和本发明中的a‘₂×b’₂，a‘₂×b′₃，a‘₃×b’₂以及a‘₃×b’₃。 不同之处在于，本发明至多需要4个减法来得到参数a‘₂，a‘₃，b’₂以及b’₃，而原始算法需要 5个乘法来得到a₁×b₂×2^-2(mod p)，a₁×b₂，a₁×b₃，b₁×a₂以及b₁×a₃。此外，在Barret除 法之前，原始算法需要6至9个加减法，2个右移操作，来计算参数c₁，c₂和c₃，而本发 明只需2个加法操作和1个右移操作。原始算法需要预计算2^-2(mod p)，而本发明不需要。最重要的是，在本发明中，在表示A′与B′时，权重为R²的这一项已经完全消失了， 这大大节省了操作数的数目。

三、本发明的硬件实现

为进一步说明本发明的技术方案及技术效果，本实施例提供了一个可以用于实现本 发明的硬件结构，如图5所示。该结构主要由N/2比特的乘法器、5N/2比特的加法器 以及2N比特的减法器构成。寄存器_A中存放了A′和B′等初始数据，寄存器_B和寄存 器_C分别存放了计算过程中的中间值以及最后结果。输入位宽的选择是为了保证尽可能 减少算法所需时钟数。整个模乘计算过程由有限状态机控制。N表示位宽长度，P为常 值输入端，即此处输入的为图3中P的值。位于加法器与后级选择器之间的操作代表的 是左移N/2比特操作。

在本发明，存在着3种输入长度的普通乘法操作，它们分别为N×N/2，N×N以 及3N/2×3N/2，而在输入长度为3N/2×3N/2的乘法中，其中一个乘数为常数，并且 高N/2位一直等于2，因此，可以将该乘法转换成一个输入长度为3N/2×N的乘法以 及移位操作。因此对于以上3种乘法操作，可以用同一种流水线结构实现，如图5所示。 根据图5所示结构可得，对于所述任意一种长度的乘法，在第2个时钟便可以得到结果 的最高N/2位；每当下一个时钟到达时，最终结果的下一个低N/2位便可依次得到。 例如，只需要5个时钟，我们便可以得到N×N的完整结果。

本实施例利用Vivado14.6在Virtex 7 FPGA(xc7k325tffg900-2)上构建了上述硬件结构，并选取p＝2·2³⁸⁶3²⁴²-1。

将本实施例的参数结果与现有技术的参数结果进行比对，比对的结果如表1所示：

表1

方法	FFs	LUTs	DSP48s	频率(MHz)	时钟(s)	总时间(μs)
							EFFM	11924	12790	0	31	236	7.61
本发明	9675	16629	122	55	64	1.16

从表中可以看出，本实施例的电路结构频率为55MHz，完成一次完整模乘所需时钟数 为64个，总时间为1.16us。在资源消耗增加不多的情况，吞吐量达到了原来的6倍多。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员 来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种SIDH中特殊域快速模乘的方法，其特征在于，包括SIDH中单个模乘的计算步骤：

(1)定义SIDH中单个模乘为A·B(mod p)，其中，p为素数域，A和B为素数域p中的任意两个元素；p＝2·2^a3^b-1，a，b均为偶数；

(2)将A·B(mod p)表示为：

式中，A、A′、B、B＇分别用以R为基的补码形式表示为：

A＝a₁R²+a₂R+a₃

B＝b₁R²+b₂R+b₃

A′＝a′₂R+a′₃

B′＝b′₂R+b′₃

其中，R＝2^a/23^b/2，a₁、a₂、a₃、b₁、b₂、b₃、a′₂、a′₃、b′₂、b′₃均为系数，符号表示异或；

(3)通过计算A′·B′(mod p)得到A·B(mod p)，包括依次执行步骤(3-1)至(3-8)：

(3-1)将A′·B′(mod p)的值以基R形式表示为：

A′·B′(mod p)＝c′₁R²+c′₂R+c′₃

式中，c＇₁、c′₂、c′₃为系数变量；

(3-2)初始化c′₁＝0，c′₂＝0，c′₃＝0；

(3-3)计算c′₁＝c′₂c′₃[0]，[0]表示取最低位；

(3-4)计算

(3-5)对c′₃进行Barrett约减，约减结果记为r₁，即c′₃＝r₁；计算c′₂＝c′₂+r₁；

(3-6)对c′₂进行Barrett约减，约减结果记为r₂，即c′₂＝r₂；计算c＇₁＝c＇₁+r₂；

(3-7)计算c′₁＝c＇₁[0]；

(3-8)根据计算得到的c′₁、c′₂、c′₃得到A′·B′(mod p)的计算结果，再根据A′·B′(modp)得到A·B(mod p)的计算结果。