clmul

contributed by < laochanlam >, < etc276 >, < jack81306 >, < peterting >

tags: laochanlam etc276

相關資訊

開發環境

• Ubuntu 16.10 ( 64 bit )

Architecture:          x86_64
CPU op-mode(s):        32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                4
On-line CPU(s) list:   0-3
Thread(s) per core:    2
Core(s) per socket:    2
Socket(s):             1
NUMA node(s):          1
Vendor ID:             GenuineIntel
CPU family:            6
Model:                 61
Model name:            Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz
Stepping:              4
CPU MHz:               2400.878
CPU max MHz:           3000.0000
CPU min MHz:           500.0000
BogoMIPS:              4788.90
Virtualization:        VT-x
L1d cache:             32K
L1i cache:             32K
L2 cache:              256K
L3 cache:              4096K
NUMA node0 CPU(s):     0-3

相關連結

• 2017 年春季作業說明
• 2017 年春季班第一次分組表
• 2017 年春季班第二次分組表
• B04: clz 作業要求
• 課程進度與開放資源 Wiki
• Youtube 解說
• Github

作業要求

延伸 B04: clz 要求，探討 carryless multiplication 演算法和實作 (需要考慮到 Intel 延伸指令集)，介紹密碼學的應用

背景知識

• 加法器與乘法器
• carry-less mutiplication (無進位乘法)
• Galois Fields
• Intel® AES NI
• GCM & AES
• PCLMULQDQ 在 ECC 的應用

---

加法器與乘法器

• 加法器是用於執行加法的電路元件，通常由 and閘、or閘 和 xor閘構成
  • 也可以用加法器實現減法，只要將減數轉成二補數，並注意溢位即可
  • 半加器：將兩個一位二進位數相加 (input: A, B) (output: S, C)
    
    
  • 全加器：將兩個一位二進位數相加 (input: A, B, Cin) (output: S, Cout)
    
    
  • 波紋進位加法器：使用多個一位全加器構成 N 位加法器
    
    
  • 半加器簡易的 code 可表示為

uint half_add(uint a, uint b)
{
    if (b == 0) return a;
    uint sum = a ^ b;             /* 相加但不進位 */
    uint carry = (a & b) << 1;    /* 進位但不相加 */
    return half_add(sum, carry);
}

• 乘法器則是利用加法器，依序判斷乘數的最後一位決定是否相加，並在每次迴圈中都左/右移
  • 用 c code 表示

uint mul(uint a, uint b)
{
    uint r = 0;
    while (b != 0)
    {
        if ((b & 1) != 0)
            r += a;                /* 執行加法 */
        a <<= 1;
        b >>= 1;
    }
    return r;
}

carry-less mutiplication (無進位乘法)

• 加法的概念中，包含相加 (xor) 和進位 (and)
• 無進位乘法類似一般的乘法，差別僅在於將 + 改為 ^ (捨去進位)
• 用 c code 表示 (跟上方的 code 僅有一行的差異)

uint mul(uint a, uint b)
{
    uint r = 0;
    while (b != 0)
    {
        if ((b & 1) != 0)
            r ^= a;                /* 執行無進位加法 */
        a <<= 1;
        b >>= 1;
    }
    return r;
}

• 以 3乘以3 為例 (x表示乘法，*表示無進位乘法)
  3 x 3 = 9
  3 * 3 = 5

    11        11
x   11    *   11
------    ------
    11        11
+  11     ^  11
------    ------
  1001       101 

Galois Fields

• 體（Field）
  • 一種可進行加、減、乘和除運算的代數結構。
  • 體有這樣一個性質：在加法和乘法上具有封閉性。
  • 但是要注意，體裡面的乘法和加法不一定是我們平常使用的乘法和加法。
  • 一個集合有加法單位元素和乘法單位元素
  • 每一個元素都對應有加法反元素和乘法反元素，是成為體的必要條件。

選用適當的數學結構可保證計算上的正確性，以及利用其數學性質將運算加速

• 有限體 (Finite Field 或 Galois Field)
  是包含有限個數元素的體。有限體最常見的例子是當 p 為質數時，整數對 p 取模。

• 以 GF(5) 為例
  GF(5) = {0, 1, 2, 3, 4}
  2 + 4 = 1 (mod 5)，亦即 2 的加法反元素是 4
  2 * 3 = 1 (mod 5)，亦即 2 的乘法反元素是 3

• p 為質數時，能保證集合中的所有的元素都有加法和乘法反元素 (0除外)

  • show :
    若 p 不為質數，則 0~p 之間必定存在一個 a，使得 a, p 不互質，而不存在 a 的乘法反元素
  • proof :
    歐拉定理證明了當 a 和 n 互質時，會滿足
    這條方程式。φ(n) 是歐拉函數，代表小於 n 且和 n 互質的個數。因此上述結果可分解為
    其中的
    即為 a 對於 n 的模反元素。若 n 不為質數時，則
    中的指數部分可能為 0 ,因此若要必定存在模反元素，則 n 必為質數。

• 另一種常見的有限體 GF(pⁿ)

  • 乘法和加法均為無進位 (+ 改為 xor)
  • 將集合內元素表達為 f(x) = a_n-1x^n-1+a_n-2x^n-2+…+a₁x+a₀, a_i∈[0,p)
  • 維持封閉性的方法為設定一個質多項式取模
  • 特別的當 p = 2 時，加法和減法等價 (+1 = -1 = xor 1)

• 以 GF(2⁸) 為例
  使用 m(x) = x⁸ + x⁴ + x³ + x² + 1 作為質多項式

取模的概念是，將到達上界的值替換掉
以 2 * 6 (mod 10) 為例，就是計算 2 * 6 = 12 之後，將 10 替換為 0 ，得到 2
所以這邊就是將所有的 x⁸ 替換成 x⁴ + x³ + x² + 1 (正負號相同)

GF(28)	a的多項式表達	D7D6D5D4D3D2D1D0	推導過程
0	0	0 0 0 0 0 0 0 0
a0	a0	0 0 0 0 0 0 0 1
a1	a1	0 0 0 0 0 0 1 0
a2	a2	0 0 0 0 0 1 0 0
a3	a3	0 0 0 0 1 0 0 0
a4	a4	0 0 0 1 0 0 0 0
a5	a5	0 0 1 0 0 0 0 0
a6	a6	0 1 0 0 0 0 0 0
a7	a7	1 0 0 0 0 0 0 0
a8	a4+a3+a2+a0	0 0 0 1 1 1 0 1	P(a) = a8 + a4 + a3 + a2 + a0
a9	a5+a4+a3+a1	0 0 1 1 1 0 1 0	a8*a = (a4+a3+a2+a0)*a = a5+a4+a3+a1
a10	a6+a5+a4+a2	0 1 1 1 0 1 0 0	a9*a = (a5+a4+a3+a1)*a = a6+a5+a4+a2
a11	a7+a6+a5+a3	1 1 1 0 1 0 0 0	a10*a = (a6+a5+a4+a2)*a = a7+a6+a5+a3
a12	a7+a6+a3+a2+a0	1 1 0 0 1 1 0 1	a11*a = (a7+a6+a5+a3)*a = a7+a6+a3+a2+a0
a13	a7+a2+a1+a0	1 0 0 0 0 1 1 1	a12*a = (a7+a6+a3+a2+a0)*a = a8+a7+a4+a3+a1 = (a4+a3+a2+a0) + a7+a4+a3+a1 = a7+a2+a1+a0

在 GF(pⁿ)的情況下，必須找一個 n 次方的質多項式 (不被集合內的元素整除)，以確保封閉性

如同上面所說，把超過上界(pⁿ)的部分進行替換，且因為沒有進位的概念，不會替換完又超過 pⁿ，就可以確保所有運算均具有封閉性

• 以數字表示

  • 以 GF (16) 為例
    若 A = 14, B = 11 則 mul (A, B) = 10
  ​​​​C = A * B = 154 = 9 * 16 + 10
  

  • 以 GF (2⁴) 為例，質多項式為 x⁴+x+1
    若 A = [1110], B = [1011]，則 cl_mul(A, B) = [1000]
  ​​​​C = A∙B = [01100010] = [0110]∙[10011] XOR [1000]
  

• 以多項式表示

  3 * 3
  = (x+1)(x+1)
  = x² + 2x + 1 = 9₍₁₀₎ (進位)
  = x² + 1 = 5₍₁₀₎ (無進位)

• 無進位乘法在密碼學上的應用

  • const time : 避免時間差攻擊 (本頁下方有實驗)
  • 雪崩效應 : xor 相對於 and 和 or，有 50% 的機率出現反轉
  • 非線性 : 運算的規則較無規律，不易反推

參考及引用資料

資訊安全 Galois Fields C++ 軟體實作體驗
伽罗华域(Galois Field，GF，有限域)乘法运算

---

Intel® 進階加密標準新增指令 (Intel® AES NI)

在尋找 Intel 有關加密的延伸指令集時，找到了 AES-NI 的指令集，而在官方的文件中看到了以下一段：

PCLMULQDQ – Performs carry-less multiplication of two 64-bit data into a 128-bit result.
Carry-less multiplication of two 128-bit data into a 256-bit result can use PCLMULQDQ as building blocks.
Carry-less multiplication is an important piece of implementing Galois Counter Mode (GCM) operation of block ciphers.
GCM operation can be used in conjunction with AES algorithms to add an authentication capability.
GCM usage models also include IPsec, storage standard, and security protocols over fiber channel. Additionally, PCLMULQDQ can be used in
calculations of hash functions and CRC using arbitrary polynomials.

有提及到該指令集中有一個指令 PCLMULQDQ ，可輸入兩個 64-bit 的 data 作 Carry-less multiplication 輸出 128-bit 的 data。

而此處亦有提及到 Carry-less muliplication 是 GCM 加密應用的重要一環，所以先來研究一下 GCM (Galois/Counter Mode)。

GCM (Galois/Counter Mode) 的介紹在這裡。

隨後找到了 PCLMULQDQ instruction 的官方文件，先來把文件讀完，隨後進行實作！

PCLMULQDQ instruction 官方文件導讀

• 簡介
  • PCLMULQDQ 是 intel 的處理器指令，可以計算 2 個 64 bit 的無進位乘法
  • 這份 paper 提供一些指令細節和各種演算法 for GCM
  • 而這可以應用在 AES-GCM 上
• 基礎知識
  • 無進位乘法是許多加密的基礎，所以加快無進位乘法很重要
  • 無進位乘法的定義
  • 無進位乘法在 Galois Field 的情況
• 內文
  • PCLMULQDQ Instruction Definition (和 AVX 版的)
  • The Galois Counter Mode (GCM，下方有介紹)
  • 執行無進位乘法的演算法 (algorithm 1, 2)
  • 有效率的 reduced algorithm (一些數學的推導)
  • 應用在 x¹²⁸ + x⁷ + x² + x + 1 的 reduce
  • Bit Reflection Peculiarity (有點像字串反轉)
  • 一些可以避免反轉的替代方案 (因為反轉不只是 endian 的問題)
  • 實作 (Using Linear Folding)
• 部分說明
  • 可以用 4 to 4，也就是 8 個 bit 記錄下 reflection (查表比較快)
  • xmm1 紀錄 low bit
  • xmm2 紀錄 high bit
  • xmm3 紀錄 low bit 的 reflec (e.g. 0->0, 1->8)
  • xmm4 紀錄 high bit 的 reflec

參考及引用資料

進階加密指令集和安全金鑰構成 Intel® 資料保護技術 (Intel® Data Protection Technology)
INTRODUCTION TO INTEL® AES-NI AND INTEL® SECURE KEY INSTRUCTIONS
CPUIDFIELD：CPUID字段的统一编号、读取方案。范例：检查SSE4A、 AES、PCLMULQDQ指令
Intel® Carry-Less Multiplication Instruction and its Usage for Computing the GCM Mode
AES指令集 Wiki

---

GCM 與 AES 的關係

關於 GCM 的詳細介紹

關於 AES 的詳細介紹

這是 GCM 加密的流程圖：

問：那麼，到底 GCM 和 AES 之間有什麼關係，以致 GCM 裡面所用到的無進位乘法要放在 AES-NI 的指令集裡呢？

• 因為 GCM 的加密流程會搭配 AES 的加密函式，合稱 AES-GCM，也就是把 AES 套用在上圖中的
  中，才能成為完整的加密流程，所以需要使用 AES 進行加密的同時，也要使用到無進位乘法來實現加密流程(GCM)。

參考及引用資料

AES-GCM 加密算法
分組密碼工作模式

---

PCLMULQDQ 在 ECC 的應用

• ECC 是一種建立公開金鑰加密的演算法，基於橢圓曲線數學的離散對數難題而成立

• 簡易的 ECC 說明

  • 選定一個橢圓曲線
  • 選定兩個線上的點，進行標量加法 (畫線找曲線上的另一個交點)
  • 若兩點相同，則是找此點的切線與曲線上的交點
  • 對稱於 x軸的交點，即為加法運算後的答案
  • 乘法為重複加法的次數，如 P+P = 2*P
    
    

• 橢圓曲線(連續) 藉由 GF 的性質而離散化

• 以 GF(5) 的 y² = x^3 + x + 1 為例

  • 曲線上面的點有(0,1) (0,4) (2,1) (2,4) (3,1) (3,4) (4,3) (4,2) 和無窮遠，共九個點

• 離散對數難題

  • P, Q 為橢圓線上的點，k為常數
  • 給定 P 和 k，可輕易求得 Q = k*P
  • 但給定 P 和 Q，難以在時間內計算出 k

• 其他細節

  • 有限體的選擇 (GF(2ⁿ)在硬體上實現較快)
  • 橢圓曲線的選擇 (攸關密碼的長度及安全性)
  • 點的表達方式 (例如 L´opez-Dahab projective coordinates)
  • 計算上相關的演算法

• Montgomery Ladder

  • 基本的乘冪為單純從左至右或從右至左相乘
  • 存在 condition jump，會被攻擊 (時間分析、電源不穩等等)
  • 概念為每次計算都為了下一次，即使有分支也做相似的事情
  • 以下是虛擬碼

  R0 ← 0
  R1 ← P
  for i from m downto 0 do
     if di = 0 then
        R1 ← point_add(R0, R1)
        R0 ← point_double(R0)
     else
        R0 ← point_add(R0, R1)
        R1 ← point_double(R1)
  return R0

• Karatsuba
  • 兩個 n 位數相乘，原本需要 O(n^2) 的時間
  • 此演算法可以降成 T(n) = 3*T(n/2) + cn +d ，也就是 O(n^ log3)
  • 可應用於硬體上，使得無進位乘法更快 (可參考 這篇 )

ECC 的優點是在相同安全度下，所需的密鑰長度較短，但缺點就是計算較慢，所以利用 PCLMULQDQ 來加速，參考資料 [2] 中提到這樣可以提升 10% 的速度

參考資料

Intel Polynomial Multiplication Instruction and its Usage for Elliptic Curve Cryptography
Software implementation of binary elliptic curves:impact of the carry-less multiplier on scalar multiplication
近代加密 快速冪次計算
橢圓曲線密碼系統

開發紀錄

開發目標

1. [x]設計實驗測試 Intel 延伸指令集的優化程度

2. [x]設計 Intel 延伸指令集的數值分佈實驗

首先這個 AES-NI 是一個 Instruction set 嘛，所以我來找找怎麼把他放到 c code 裡面。

<mmintrin.h>  MMX

<xmmintrin.h> SSE

<emmintrin.h> SSE2

<pmmintrin.h> SSE3

<tmmintrin.h> SSSE3

<smmintrin.h> SSE4.1

<nmmintrin.h> SSE4.2

<ammintrin.h> SSE4A

<wmmintrin.h> AES

<immintrin.h> AVX

<zmmintrin.h> AVX512

這個表格不錯，先把他記下來。

後來在參考完這裡之後，得知若要用 C 寫 Instruction Set 有以下兩種方法：

1. 內嵌式組合語言
   就是很醜的把 Instruction 內嵌進去程式碼裡面，像這樣：

    __asm addps xmm0, xmm1
    __asm movaps [ebx], xmm0
    ...
    __m128 data;
    ...
    __asm
    {
        lea ebx, data
        addps xmm0, xmm1
        movaps [ebx], xmm0
    }

2. 使用 intrinsics
   就像我們之前使用 SSE AVX 等的方法如下：

    __m128 data1, data2;
    ...
    __m128 out = _mm_add_ps(data1, data2);
    ...

雖然後者的可讀性比較高，可是還是要經過 Compiler 的產生，不同的 Compiler 有可能會對效能有不同的影響。

接下來要開始去找 AES-NI 的 intrinsics 了！

先是從 Intel® Advanced Encryption Standard Instructions (AES-NI) 看到這一段

Using C/C++ or assembly

If you have existing C/C++ or assembly implementations of AES algorithms you can take advantage of the support 
provided in most of the standard compiler development tools. You will need to modify your code to replace code 
blocks with the equivalent AES-NI instructions. AES-NI instructions can be called from C/C++ either using 
inline assembly or using special functions know as intrinsics. Each intrinsic maps to one of the new 
instructions. Using intrinsics allows you to develop code using the syntax of C/C++ function calls and 
variables instead of inline assembly language.

To use AES-NI in assembly language you can directly call the relevant instruction from your code.

The following compilers provide C/C++ as well assembly support for AES-NI.

Compiler	Description	Version supporting AES-NI
Gcc/g++	Open source GNU compiler for C/C++	4.4 or later
Intel® C/C++ compiler	Intel compiler tools for C/C+	11.1 or later
Microsoft* Visual C++	C/C++ compiler tools for Windows* operating systems	2008 SP1 or later

意思呢，大概就是說可以用 intrinsics 實現 AES-NI 的 Instruction Set。

之後就找到了 wmmintrin.h ( intrinsics of AES-NI ) 的 Document，來看看 Source Code：

    /*===---- wmmintrin.h - AES intrinsics ------------------------------------===
      *
      * Permission is hereby granted, free of charge, to any person obtaining a copy
      * of this software and associated documentation files (the "Software"), to deal
      * in the Software without restriction, including without limitation the rights
      * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
      * copies of the Software, and to permit persons to whom the Software is
      * furnished to do so, subject to the following conditions:
      *
      * The above copyright notice and this permission notice shall be included in
      * all copies or substantial portions of the Software.
      *
      * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
      * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
      * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
      * THE SOFTWARE.
      *
      *===-----------------------------------------------------------------------===
      */
      
     
     #ifndef _WMMINTRIN_H
     #define _WMMINTRIN_H
     
     #include <emmintrin.h>
     
     #include <__wmmintrin_aes.h>
     
     #include <__wmmintrin_pclmul.h>
     
     #endif /* _WMMINTRIN_H */

有我們一直在找的 clmul ，然後依照他給的的文件進行實作。

指令實作

• 讓我們先來看 _mm_clmulepi64_si128 這個指令
  • macro def : __m128i _mm_clmulepi64_si128(__m128i __X, __m128i __Y, const int __I)
  • X 和 Y 都是 __mm128i (由兩個 64bit vector 所組成的) parameter
  • const 則是代表我們的 clmul 輸出是由哪個部分的 X, Y 所組成
  • 而完整的 clmul (X, Y) 會是 256 bit
  • 但我們只需要 128 bit 所以 X, Y 各只會挑一段來用
  • 值得一提的是這個指令對於第三個參數 only looks at the least significant bit of each

example code :

#include <wmmintrin.h>
#include <stdio.h>

int main()
{
    __m128i a;
    __m128i b;

    a.m128i_i64[1] = 2;
    a.m128i_i64[0] = -1284;
    b.m128i_i64[1] = 25;
    b.m128i_i64[0] = 65535;

    // _mm_clmulepi64_si128 only looks at the least significant bit of each
    // hexadecimal integer
    const int product1 = 0x11;
    const int product2 = 0x00;
    const int product3 = 0xF2;

    int expect1 = int ( a.m128i_i64[1] * b.m128i_i64[1] );
    int expect2 = int ( a.m128i_i64[0] * b.m128i_i64[0] );
    int expect3 = int ( a.m128i_i64[0] * b.m128i_i64[1] );

    __m128i result1 = _mm_clmulepi64_si128( a, b, product1 );
    __m128i result2 = _mm_clmulepi64_si128( a, b, product2 );
    __m128i result3 = _mm_clmulepi64_si128( a, b, product3 );

    printf_s("%I64d times %I64d without a carry bit: %I64d\n", 
        a.m128i_i64[1], b.m128i_i64[1], result1.m128i_i64[0]);
    printf_s("%I64d times %I64d without a carry bit: %I64d\n", 
        a.m128i_i64[0], b.m128i_i64[0], result2.m128i_i64[0]);
    printf_s("%I64d times %I64d without a carry bit: %I64d\n", 
        a.m128i_i64[0], b.m128i_i64[1], result3.m128i_i64[0]);

    return 0;
}

錯誤訊息紀錄 :

2 times 25 without a carry bit: 50
-1284 times 65535 without a carry bit: 50419284
-1284 times 25 without a carry bit: -32036

在執行時我們遇到一個問題，將 a, b 宣告成 __m128i 但是卻得到錯誤訊息，

request for member ‘m128i_i64’ in ‘a’, which is of non-class type ‘__m128i {aka __vector(2) long long int

然後發現我不小心把 .c 檔打成 .cpp (眼神死)。

改回 .c 檔出現以下錯誤訊息，

error: request for member ‘m128i_i64’ in something not a structure or union

查了一下發現 gcc 會直接把 __m128i 當成 long[2]，所以不支援.m128i_i64[1]這類的 function ，所以手動改成 a[1] = 2 。

又得到以下錯誤，

error: ‘__builtin_ia32_pclmulqdq128’ needs isa option -m32 -mpclmul

補上編譯選項之後，

fatal error: sys/cdefs.h: No such file or directory #  include <sys/cdefs.h>

砍掉 -m32 的編譯選項後，

error: the last argument must be an 8-bit immediate

將 product 從 int 改成 int8_t 仍不能解決問題, 最後將第三個參數直接放 0x00 才可以跑。

以下是最後可以跑的版本

#include <wmmintrin.h>
#include <stdio.h>

int main()
{
    __m128i a;
    __m128i b;

    a[1] = 2;
    a[0] = -1284;
    b[1] = 25;
    b[0] = 65535;

    // _mm_clmulepi64_si128 only looks at the least significant bit of each

    int expect1 = a[1] * b[1];
    int expect2 = a[0] * b[0];
    int expect3 = a[0] * b[1];

    __m128i result1 = _mm_clmulepi64_si128( a, b, 0x11 );
    __m128i result2 = _mm_clmulepi64_si128( a, b, 0x00 );
    __m128i result3 = _mm_clmulepi64_si128( a, b, 0xF2 );

    printf("%I64d times %I64d without a carry bit: %I64d\n",
        a[1], b[1], result1[0]);
    printf("%I64d times %I64d without a carry bit: %I64d\n",
        a[0], b[0], result2[0]);
    printf("%I64d times %I64d without a carry bit: %I64d\n",
        a[0], b[1], result3[0]);
    return 0;
}

---

優化程度實驗

為了得知 _mm_clmulepi64_si128 在效能上的提升，我們設計了以下實驗：

1. 使用 Intrinsics 的版本：

        __m128i I0 = _mm_loadu_si128((__m128i *)(src));

        __m128i result1 = _mm_clmulepi64_si128( I0, I0, 0x00 ); // [63:0] [63:0]
        __m128i result3 = _mm_clmulepi64_si128( I0, I0, 0xF2 ); // [63:0] [127:64]
        __m128i result2 = _mm_clmulepi64_si128( I0, I0, 0x01 ); // [127:64] [63:0]
        __m128i result4 = _mm_clmulepi64_si128( I0, I0, 0xFF ); // [127:64] [127:64]

        _mm_storeu_si128((__m128i *)(dst1), result1);
        _mm_storeu_si128((__m128i *)(dst2), result2);
        _mm_storeu_si128((__m128i *)(dst3), result3);
        _mm_storeu_si128((__m128i *)(dst4), result4);

2. 經修改後 cl_mul 實作的版本：

static int64_t cl_mul(int64_t a, int64_t b)
{
    int64_t r = 0;
    while (b != 0)
    {
        if ((a & 1) != 0)
            r ^= b;      // carryless addition is xor
        a >>= 1;
        b <<= 1;
    }
    return r;
}

我們把以上兩個版本分別進行 100 次的運算，計算出來平均值的繪圖如下：

可見有使用 Instruction Set ( AES-NI ) 實作的 clmul 比原來的程式效能提升有超過十倍之多。

實驗準確性

在實驗的準確性方面，我們使用了 assert 來進行驗證。

    assert( cl_mul(random_number[0],random_number[0]) == dst1[0]);
    assert( cl_mul(random_number[1],random_number[0]) == dst2[0]);
    assert( cl_mul(random_number[0],random_number[1]) == dst3[0]);
    assert( cl_mul(random_number[1],random_number[1]) == dst4[0]);

我們對 _mm_clmulepi64_si128 的四種情況各進行了驗證，以確定實驗得出的答案一致。

可是我們遇到了一個問題，就是 _mm_clmulepi64_si128 會返回 128-bit，但我們實作的 cl_mul 只能做到返回 64-bit，若改為用數組作儲存方式的話時間會大大增加，所以我們只有作 _mm_clmulepi64_si128 的後 64-bit 直接與 cl_mul 的 64-bit 作比對。
laochanlam

---

數值分佈實驗

接下來為了實驗指令 _mm_clmulepi64_si128(a, b, product) 是否在效能上達到 const time ，設計了以下實驗：

因為無進位乘法需要的參數有 a, b 兩個，所以以下分成幾個狀況
laochanlam

a, b 為 10000 個隨機數字分佈:

程式碼:

    for (int j = 0; j < 100000; j++) {
        for (int i = 0; i < 2; i++)
            src[i] = ( rand() % 922337203685477580 ) + 1;
        clock_gettime(CLOCK_REALTIME, &start);

        for (int k = 0; k < 10000; k++) {
            __m128i I0 = _mm_loadu_si128((__m128i *)(src));

            __m128i result1 = _mm_clmulepi64_si128( I0, I0, 0x00 ); // [63:0] [63:0]
            __m128i result3 = _mm_clmulepi64_si128( I0, I0, 0xF2 ); // [63:0] [127:64]
            __m128i result2 = _mm_clmulepi64_si128( I0, I0, 0x01 ); // [127:64] [63:0]
            __m128i result4 = _mm_clmulepi64_si128( I0, I0, 0xFF ); // [127:64] [127:64]

            _mm_storeu_si128((__m128i *)(dst1), result1);
            _mm_storeu_si128((__m128i *)(dst2), result2);
            _mm_storeu_si128((__m128i *)(dst3), result3);
            _mm_storeu_si128((__m128i *)(dst4), result4);
        }

        clock_gettime(CLOCK_REALTIME, &end);
    }

a = b 時遞增的數值分佈效能狀況:

    for (int j = 0; j < 100000; j++) {
            src[0] = value;
        clock_gettime(CLOCK_REALTIME, &start);

        for (int k = 0; k < 10000; k++) {
            __m128i I0 = _mm_loadu_si128((__m128i *)(src));

            __m128i result4 = _mm_clmulepi64_si128( I0, I0, 0x00 ); // [63:0] [63:0]
            _mm_storeu_si128((__m128i *)(dst4), result4);
        }
        value += 10000;
        clock_gettime(CLOCK_REALTIME, &end);
    }

當 a 設定為常數，分析在 b 在數值改變下的效能狀況：

    src[1] = 1234567896;
    for (int j = 0; j < 100000; j++) {
            src[0] = value;
        clock_gettime(CLOCK_REALTIME, &start);

        for (int k = 0; k < 10000; k++) {
            __m128i I0 = _mm_loadu_si128((__m128i *)(src));

            __m128i result3 = _mm_clmulepi64_si128( I0, I0, 0xF2 ); // [63:0] [127:64]
            _mm_storeu_si128((__m128i *)(dst3), result3);

        }
        value += 10000;
        clock_gettime(CLOCK_REALTIME, &end);
    }

由以上實驗可得出結論 _mm_clmulepi64_si128 的指令為 Const Time。

參考及引用資料整理：

• ECC

  • Intel Polynomial Multiplication Instruction and its Usage for Elliptic Curve Cryptography
  • Software implementation of binary elliptic curves:impact of the carry-less multiplier on scalar multiplication
  • 近代加密 快速冪次計算
  • 橢圓曲線密碼系統

• Galois Fields

  • 資訊安全 Galois Fields C++ 軟體實作體驗
  • 伽罗华域(Galois Field，GF，有限域)乘法运算

• carryless multiplication

  • carryless multiplication
  • Carry-Less Multiplication Instruction, Usage for the GCM Mode

• AES-NI

  • 進階加密指令集和安全金鑰構成 Intel® 資料保護技術 (Intel® Data Protection Technology)
  • INTRODUCTION TO INTEL® AES-NI AND INTEL® SECURE KEY INSTRUCTIONS
  • CPUIDFIELD：CPUID字段的统一编号、读取方案。范例：检查SSE4A、 AES、PCLMULQDQ指令
  • Intel® Carry-Less Multiplication Instruction and its Usage for Computing the GCM Mode
  • AES指令集 Wiki

• AES-NI

  • AES-GCM 加密算法
  • 分組密碼工作模式

• _mm_clmulepi64_si128 指令

  • Intel® Advanced Encryption Standard Instructions (AES-NI)
  • wmmintrin.h File Reference
  • _mm_clmulepi64_si128
  • wmmintrin.h
  • SSE 介紹