384位NIST素域椭圆曲线快速约减算法x64编程实现研究(上)

最新推荐文章于 2023-08-02 20:20:43 发布

weixin_34185560

最新推荐文章于 2023-08-02 20:20:43 发布

阅读量671

点赞数

文章标签： c/c++ python

原文链接：https://my.oschina.net/safedead/blog/377237

版权

2019独角兽企业重金招聘Python工程师标准>>>

384位NIST素域椭圆曲线的公开资料中对于快速约减算法是这么描述的：

摘自：Mathematical routines for the NIST prime elliptic curves (April 05, 2010)

p384 = 2 ^ 384 − 2 ^ 128 − 2 ^ 96 + 2 ^ 32 − 1
p384 = ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff
       ffffffff fffffffe ffffffff 00000000 00000000 ffffffff

Routine 3.2.11 mp_mod_384 (r, a): Set r = a (mod p384 )
1: {Note: the ai are 32–bit quantities.}
2:  t  ← ( a11 || a10 || a9  || a8  || a7  || a6  || a5  || a4  || a3  || a2  || a1  || a0  )
3:  s1 ← (  0  ||  0  || 0   || 0   || 0   || a23 || a22 || a21 ||  0  ||  0  ||  0  ||  0  )
4:  s2 ← ( a23 || a22 || a21 || a20 || a19 || a18 || a17 || a16 || a15 || a14 || a13 || a12 )
5:  s3 ← ( a20 || a19 || a18 || a17 || a16 || a15 || a14 || a13 || a12 || a23 || a22 || a21 )
6:  s4 ← ( a19 || a18 || a17 || a16 || a15 || a14 || a13 || a12 || a20 ||  0  || a23 ||  0  )
7:  s5 ← (  0  ||  0  ||  0  ||  0  || a23 || a22 || a21 || a20 ||  0  ||  0  ||  0  ||  0  )
8:  s6 ← (  0  ||  0  ||  0  ||  0  ||  0  ||  0  || a23 || a22 || a21 ||  0  ||  0  || a20 )
9:  d1 ← ( a22 || a21 || a20 || a19 || a18 || a17 || a16 || a15 || a14 || a13 || a12 || a23 )
10: d2 ← (  0  ||  0  ||  0  ||  0  ||  0  ||  0  ||  0  || a23 || a22 || a21 || a20 ||  0  )
11: d3 ← (  0  ||  0  ||  0  ||  0  ||  0  ||  0  ||  0  || a23 || a23 ||  0  ||  0  ||  0  )
12: d1 ← p384 − d1
13: r ← t + 2s1 + s2 + s3 + s4 + s5 + s6 + d1 − d2 − d3
14: Reduce r mod p384 by subtraction of up to four multiples of p384.

为了用x64汇编语言实现这个算法，需要对算法的描述形式加以转变，使之更适用于编程。

第1阶段转换：水平翻转数值矩阵使其变为先低位后高位的表达形式

|a00|a01|a02|a03|a04|a05|a06|a07|a08|a09|a10|a11| ->t
| 0 | 0 | 0 | 0 |a21|a22|a23| 0 | 0 | 0 | 0 | 0 | ->s1
| 0 |a23| 0 |a20|a12|a13|a14|a15|a16|a17|a18|a19| ->s2
|a21|a22|a23|a12|a13|a14|a15|a16|a17|a18|a19|a20| ->s3
|a12|a13|a14|a15|a16|a17|a18|a19|a20|a21|a22|a23| ->s4
| 0 | 0 | 0 | 0 |a20|a21|a22|a23| 0 | 0 | 0 | 0 | ->s5
|a20| 0 | 0 |a21|a22|a23| 0 | 0 | 0 | 0 | 0 | 0 | ->s6
-------------------------------------------------
|a23|a12|a13|a14|a15|a16|a17|a18|a19|a20|a21|a22| ->d1
| 0 |a20|a21|a22|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 | ->d2
| 0 | 0 | 0 |a23|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 | ->d3

第2阶段转换：按照元素连续性分组

分组1:(a00 ~ a19)
|a00|a01|a02|a03|a04|a05|a06|a07|a08|a09|a10|a11|
| 0 | 0 | 0 | 0 |a12|a13|a14|a15|a16|a17|a18|a19|
| 0 | 0 | 0 |a12|a13|a14|a15|a16|a17|a18|a19| 0 |
|a12|a13|a14|a15|a16|a17|a18|a19| 0 | 0 | 0 | 0 |
-------------------------------------------------
| 0 |a12|a13|a14|a15|a16|a17|a18|a19| 0 | 0 | 0 |

分组2:(a20 ~ a23)
| 0 | 0 | 0 | 0 |a21|a22|a23| 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 |a21|a22|a23| 0 | 0 | 0 | 0 | 0 |
| 0 |a23| 0 |a20| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
|a21|a22|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |a20|
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |a20|a21|a22|a23|
| 0 | 0 | 0 | 0 |a20|a21|a22|a23| 0 | 0 | 0 | 0 |
|a20| 0 | 0 |a21|a22|a23| 0 | 0 | 0 | 0 | 0 | 0 |
-------------------------------------------------
|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |a20|a21|a22|
| 0 |a20|a21|a22|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 |a23|a23| 0 | 0 | 0 | 0 | 0 | 0 | 0 |

第3阶段转换：在保证每组运算非负前提下按64位补齐调整

分组1:(完成)
|a00|a01|a02|a03|a04|a05|a06|a07|a08|a09|a10|a11|
|a12|a13|a14|a15|a16|a17|a18|a19|a20|a21|a22|a23|
| 0 | 0 | 0 | 0 |a12|a13|a14|a15|a16|a17|a18|a19|
| 0 | 0 |a23|a12|a13|a14|a15|a16|a17|a18|a19|a20|
-------------------------------------------------
|a23|a12|a13|a14|a15|a16|a17|a18|a19|a20| 0 | 0 |

分组2:(完成)
|                     p384                      |
-------------------------------------------------
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |a21|a22|

分组3:(等待继续优化)
| 0 | 0 | 0 | 0 |a21|a22|a23| 0 |
| 0 | 0 | 0 | 0 |a21|a22|a23| 0 |
| 0 |a23| 0 |a20| 0 | 0 | 0 | 0 |
|a21|a22| 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 |a20|a21|a22|a23|
|a20| 0 | 0 |a21|a22|a23| 0 | 0 |
---------------------------------
| 0 |a20|a21|a22|a23| 0 | 0 | 0 |
| 0 | 0 | 0 |a23|a23| 0 | 0 | 0 |

第4阶段转换：优化第3阶段的分组3，得到分组4和分组5

分组4:
| 0 | 0 | 0 | 0 |a20|a21|a22|a23|

分组5:
|a21|a22| 0 |a20|a21|a22|a23<<1 |
| 0 | 0 | 0 | 0 |a21|a22| 0 | 0 |
|a20|a23|a20|a21|a22|a23| 0 | 0 |
---------------------------------
| 0 |a20|a21|a22|a23<<1 | 0 | 0 |
| 0 | 0 |a20|a23| 0 | 0 | 0 | 0 |

运算顺序是先计算低位和短位长分组，计算顺序如下：

|a21|a22| 0 |a20|a21|a22|a23<<1 |(=)
| 0 |a20|a21|a22|a23<<1 | 0 | 0 |(-)
| 0 | 0 |a20|a23| 0 | 0 | 0 | 0 |(-)
| 0 | 0 | 0 | 0 |a21|a22| 0 | 0 |(+)
|a20|a23|a20|a21|a22|a23| 0 | 0 |(+)
| 0 | 0 | 0 | 0 |a20|a21|a22|a23|a20|a21|a22|a23|(+)
|a12|a13|a14|a15|a16|a17|a18|a19| 0 | 0 | 0 | 0 |(+)
| 0 | 0 | 0 | 0 |a12|a13|a14|a15|a16|a17|a18|a19|(+)
|a00|a01|a02|a03|a04|a05|a06|a07|a08|a09|a10|a11|(+)
| 0 | 0 |a23|a12|a13|a14|a15|a16|a17|a18|a19|a20|(+)
|a23|a12|a13|a14|a15|a16|a17|a18|a19|a20| 0 | 0 |(-)
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |a21|a22|(-)

运算顺序表每行末端带括号的运算符含义分别是：(=)赋值操作，即初始值；(-)本行数值符号为负，做减法；(+)表示本行符号为正，做加法。

整个流程中，唯一可能出现负数的就是最后一行做减法时，之前的累加值必然非负，可以根据各个元素出现位置证明。令a21 = a22 = (2 ^ 32 - 1)，其余元素皆为0，则累加结果得到最小值为-(2 ^ 384 - 2 ^ 320)，且此数值绝对值小于p384，根据有限域运算规则，用p384减去即可得到最终结果。若最终值大于p384，则需要减去最多不超过4个p384。

寄存器规划

数据输入指针：通用寄存器rsi

数据输出指针：通用寄存器rdi

通用寄存器备份与恢复：
|---------------|---------------|
|     xmm14     |     xmm15     |
|-------|-------|-------|-------|
|  r12  |  r13  |  r14  |  r15  |
|-------|-------|-------|-------|

数据输入缓存：
|---------------|---------------|---------------|---------------|---------------|---------------|
|     xmm0      |     xmm1      |     xmm2      |     xmm3      |     xmm4      |     xmm5      |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|a00|a01|a02|a03|a04|a05|a06|a07|a08|a09|a10|a11|a12|a13|a14|a15|a16|a17|a18|a19|a20|a21|a22|a23|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

运算过程数据：7个64位通用寄存器，从地到高排列
|-------|-------|-------|-------|-------|-------|-------|
|  r10  |  r11  |  r12  |  r13  |  r14  |  r15  |  rdx  |
|-------|-------|-------|-------|-------|-------|-------|

临时数据：5个通用寄存器
|-------|-------|-------|-------|-------|
|  rax  |  rsi  |  r8   |  r9   |  rcx  |
|-------|-------|-------|-------|-------|

其它寄存器：用于数据交换和输出到内存前的缓冲
|-------|-------|-------|
| xmm6  | xmm7  | xmm8  |
|-------|-------|-------|

简单性能测试：

1、自测程序test1.c，通过调用开源P运算库GMP的函数得到参考数据：

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#include <unistd.h>

#include <gmp.h>

int main(int argc, char *argv[])
{
        mpz_t   a, b, c, p, r;
        //int   i = 0;

        //初始化内存
        mpz_init2(a, 384);
        mpz_init2(b, 384);
        mpz_init2(c, 768);
        mpz_init2(p, 384);
        mpz_init2(r, 384);

        //初始化数据
        mpz_init_set_str(a, "EEEEEEEEEEEEEEEECCCCCCCCCCCCCCCCCD8964545891EBC5F8F0D7E9F7D5C30A3E7EB0B141E265DCC459138BCE6E7F2D", 16);
        mpz_init_set_str(b, "CD8964545891EBC5F8F0D7E9F7D5C30A3E7EB0B141E265DCC459138BCE6E7F2DEEEEEEEEEEEEEEEECCCCCCCCCCCCCCCC", 16);
        mpz_init_set_str(p, "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffeffffffff0000000000000000ffffffff", 16);

        //乘法运算
        mpz_mul(c, a, b);//c = a *b
        //求模运算
        mpz_mod(r, c, p);//r = c mod p

        //结果输出
        mpz_out_str(stdout, 16, c);
        fprintf(stdout, "\r\n");
        mpz_out_str(stdout, 16, r);
        fprintf(stdout, "\r\n");

        //释放内存
        mpz_clear(a);
        mpz_clear(b);
        mpz_clear(c);
        mpz_clear(p);

        exit(EXIT_SUCCESS);
}

编译并运行此程序，得到数据样本(过长的输出已折行)：

gcc -Wall -O2 test1.c -lgmp -o test1
./test1
bfd590d741994274668a339bec91ebef2acb6a25db704d7c7fe4faabc9b8246b
687462a249cfed730f66a066ce6f3c5d3b1a0c40603abeb865aaafee147c0410
c086d7bc4c82ac9e655392a75eef149809285ef9273c26162fb8bd29c14133dc

d0b3b0b8418585333af83fc31fdfa4f201e40c4c1b96b3a81e49401c199d271147e95cc5baf0d05e858d088f22f6feed

编写自测程序test2，引用汇编程序文件mp_mod_384.s和mul384.s，测试计算结果以及性能

#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#include <unistd.h>

#include "mul384.h"
#include "mp_mod_384.h"

int main(int argc, char *argv[])
{
        uint64_t        a[6], b[6], c[12], r[6], i;

        //初始化数据
        a[5] = 0xEEEEEEEEEEEEEEEE;
        a[4] = 0xCCCCCCCCCCCCCCCC;
        a[3] = 0xCD8964545891EBC5;
        a[2] = 0xF8F0D7E9F7D5C30A;
        a[1] = 0x3E7EB0B141E265DC;
        a[0] = 0xC459138BCE6E7F2D;

        b[5] = 0xCD8964545891EBC5;
        b[4] = 0xF8F0D7E9F7D5C30A;
        b[3] = 0x3E7EB0B141E265DC;
        b[2] = 0xC459138BCE6E7F2D;
        b[1] = 0xEEEEEEEEEEEEEEEE;
        b[0] = 0xCCCCCCCCCCCCCCCC;

        //乘法运算
        mul384(c, a, b);
        //进行1亿次求模运算
        for(i = 0; i < 100000000; i++) mp_mod_384(r, c);
        //mp_mod_384(r, c);

        //乘法结果结果输出
        for(i = 0; i < 12; i++) printf("%016lx", c[11 - i]);
        printf("\r\n");
        //求模结果输出
        for(i = 0; i < 6; i++) printf("%016lx", r[5 - i]);
        printf("\r\n");

        exit(EXIT_SUCCESS);
}

编译自测程序test2.c并计时运行：

gcc -Wall -O2 tes2.c mul384.s mp_mod_384.s -o test2
time ./test2
bfd590d741994274668a339bec91ebef2acb6a25db704d7c7fe4faabc9b8246b
687462a249cfed730f66a066ce6f3c5d3b1a0c40603abeb865aaafee147c0410
c086d7bc4c82ac9e655392a75eef149809285ef9273c26162fb8bd29c14133dc
d0b3b0b8418585333af83fc31fdfa4f201e40c4c1b96b3a81e49401c199d271147e95cc5baf0d05e858d088f22f6feed

real    0m4.395s
user    0m4.391s
sys     0m0.002s

求模结果正确。

完整汇编程序代码见后续文章：384位NIST素域椭圆曲线快速约减算法x64编程实现研究(下)

转载于:https://my.oschina.net/safedead/blog/377237