【Applied Algebra】c++大整数计算GMP库的安装和使用(Ubuntu/Mint)

c++大整数计算GMP库的安装和使用(Ubuntu/Mint)

---

GMP库(详情见在线文档)一个免费的任意精度的大数算术库,可对有符号整数、有理数和浮点数进行计算操作;除了GMP运行的计算平台中可用内存所指定的限制外,精度没有实际限制.GMP有一套丰富的功能,主要目标应用是密码学应用和研究、互联网安全应用、代数系统、计算代数研究等.

---

安装

注意加上"–enable-cxx"以使得c++可以调用(这样才会生成"libgmpxx.so");

[hanss@Mint]$ ./configure --enable-cxx
[hanss@Mint]$ make
[hanss@Mint]$ make check
[hanss@Mint]$ make install

---

编译使用

编写计算大数乘法的程序:

#include <iostream>
#include <gmpxx.h>
#include <time.h>
#include <cstdlib>
#include <ctime>
using namespace std;

int main()
{
    clock_t time = clock();
    gmp_randstate_t grt;
    gmp_randinit_default(grt);
    gmp_randseed_ui(grt, time);
    mpz_class p256(0);
    p256 = "115792089210356248762697446949407573530086143415290314195533631308867097853951";

    gmp_randclass r(gmp_randinit_default);
    mpz_class a256, b256, c256, d256;
    a256 = r.get_z_range(p256); b256 = r.get_z_range(p256); c256 = r.get_z_range(p256); d256 = r.get_z_range(p256);

    mpz_class a256_1, b256_1, c256_1, d256_1, e256_1, exp256_f, exp256_fs, exp256_leg, exp256_i;
    exp256_f = (3 * p256 - 5) / 4; exp256_fs = (p256 + 1) / 4; exp256_leg = (p256 - 1) / 2; exp256_i = p256 - 2;
    a256_1 = a256 * c256 % p256; b256_1 = b256 * d256 % p256; c256_1 = a256 * a256_1 % p256; d256_1 = b256 * b256_1 % p256; e256_1 = c256_1 + d256_1;
    mpz_t inv_1;
    mpz_init(inv_1);
    mpz_invert(inv_1, e256_1.get_mpz_t(), p256.get_mpz_t());
    mpz_class inv_2(inv_1);
    e256_1 = c256_1 * d256_1 % p256;
    e256_1 = e256_1 * inv_2 % p256;

    cout<< "exp256_f = " <<exp256_f<<" ;"<<endl;
    cout<< "exp256_fs = " <<exp256_fs<<" ;"<<endl;
    cout<< "exp256_leg = " <<exp256_leg<<" ;"<<endl;
    cout<< "exp256_i = " <<exp256_i<<" ;"<<endl;
    cout<< "e256_1 = " <<e256_1<<" ;"<<endl;
    return 0;
}

编译运行(ldconfig是为了刷新刚刚加入的动态链接库… …):

[hanss@Mint]$ g++ test.cpp -o test -lgmp -lgmpxx
[hanss@Mint]$ ldconfig
[hanss@Mint]$ ./test
exp256_f = 86844066907767186572023085212055680147564607561467735646650223481650323390462 ;
exp256_fs = 28948022302589062190674361737351893382521535853822578548883407827216774463488 ;
exp256_leg = 57896044605178124381348723474703786765043071707645157097766815654433548926975 ;
exp256_i = 115792089210356248762697446949407573530086143415290314195533631308867097853949 ;
e256_1 = 46884608919093535037063903155186519740898383003199943611564620025340708843475 ;

---

嵌入汇编的抗侧信道攻击的椭圆曲线计算

我们可以构造一个去除了"if-else"分支的cmov(嵌入汇编实现),然后再实现一系列的计算:

#include <iostream>
#include <gmpxx.h>
#include <time.h>
#include <cstdlib>
#include <ctime>
using namespace std;

inline uint64_t cmov(
    uint64_t val1,
    uint64_t test,
    uint64_t val2)
{
    uint64_t result;
    __asm__ volatile("mov %1, %0;"   // move t_val -> result
                     "test %2, %2;"  // test value of pred?
                     "cmove %3, %0;" // cmov on the test value?
                     : "=r"(result)
                     : "r"(val1), "r"(test), "r"(val2)
                     : "cc");
    return result;
}

// size: number of qwords (8 bytes)
inline void conditional_mov(
    void *source,
    void *dest,
    uint64_t size,
    uint64_t cond)
{
    __asm__ __volatile__(
        "movq %1, %%rax\n"
        "movq %2, %%rbx\n"
        "movq $0, %%rdx\n"
        "loop:\n"
        "cmpq %3, %%rdx\n"
        "jge exit\n"
        "movq (%%rbx), %%rcx\n"
        "cmp $1, %0\n"
        "cmove (%%rax), %%rcx\n"
        "movq %%rcx, (%%rbx)\n"
        "addq $8, %%rax\n"
        "addq $8, %%rbx\n"
        "inc %%rdx\n"
        "jmp loop\n"
        "exit:\n" ::"r"(cond),
        "r"(source), "r"(dest), "r"(size)
        : "rax", "rbx", "rcx", "rdx", "memory");
}

inline void __mpz_struct_cmov(__mpz_struct *a, __mpz_struct *b, uint64_t cond)
{
    conditional_mov(a->_mp_d, b->_mp_d, (a->_mp_alloc>>3)+1, cond); // 以8为单位长度计数a->_mp_alloc,取上界;
    // TODO:
    b->_mp_alloc = cmov(a->_mp_alloc, cond, b->_mp_alloc);
    b->_mp_size = cmov(a->_mp_size, cond, b->_mp_size);
}

void mpz_class_cmov(mpz_class *a, mpz_class *b, uint64_t cond) // cond 为真:把a的值给b;
{
    __mpz_struct *ma = a->get_mpz_t();
    __mpz_struct *mb = b->get_mpz_t();
    __mpz_struct_cmov(ma, mb, cond);
}

void Hush_for_Huff_f(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, //our hash encoding
                     mpz_class u, mpz_class p, mpz_class exp)
{ // f_s fast implementation, not constant time
    mpz_class v, s1, s2, gs1, gs2, x1, x2, y1, y2;
    mpz_t t;
    mpz_init(t);
    v = -u * u % p;
    s2 = -e1 * (1 + v) % p;
    gs2 = s2 * (s2 + c1) * (s2 + d1) % p;
    mpz_powm(t, gs2.get_mpz_t(), exp.get_mpz_t(), p.get_mpz_t());
    mpz_class t1(t), t2;
    t2 = t1 * t1 * gs2 % p;
    mpz_class bt, at;
    bt = b1 * t1 % p;
    at = a1 * t1 % p;
    x1 = (s2 + c1) * bt % p;
    y1 = (s2 + d1) * at % p;
 
    result[1] = u * (s2 + c1 * v % p) * bt % p;
    result[2] = u * (s2 + d1 * v % p) * at % p;

    uint64_t t3, j1, j2;
    t3 = t2.get_ui();
    mpz_class_cmov(&x1, &result[1], t3); // 抗侧信道攻击的cmov(因为消除了if-else逻辑分支);
    mpz_class_cmov(&y1, &result[2], t3);
    j1 = result[1].get_ui();
    j1 = j1 & 1;
    j2 = (j1 == t3);
    x1 = p - result[1];
    y1 = p - result[2];
    mpz_class_cmov(&x1, &result[1], j2);
    mpz_class_cmov(&y1, &result[2], j2);
}

int main()
{
    clock_t time = clock();
    gmp_randstate_t grt;
    gmp_randinit_default(grt);
    gmp_randseed_ui(grt, time);

    mpz_class p256(0);
    p256 = "115792089210356248762697446949407573530086143415290314195533631308867097853951";
    gmp_randclass r(gmp_randinit_default);

    mpz_class a256, b256, c256, d256;
    a256 = r.get_z_range(p256); b256 = r.get_z_range(p256); c256 = r.get_z_range(p256); d256 = r.get_z_range(p256);

    mpz_class a256_1, b256_1, c256_1, d256_1, e256_1, exp256_f, exp256_fs, exp256_leg, exp256_i;
    a256_1 = a256 * c256 % p256; b256_1 = b256 * d256 % p256; c256_1 = a256 * a256_1 % p256; d256_1 = b256 * b256_1 % p256; e256_1 = c256_1 + d256_1;
    mpz_t inv_1;
    mpz_init(inv_1);
    mpz_invert(inv_1, e256_1.get_mpz_t(), p256.get_mpz_t());
    mpz_class inv_2(inv_1);
    e256_1 = c256_1 * d256_1 % p256; e256_1 = e256_1 * inv_2 % p256;
    mpz_class result[2];

    // ------------ implementation for Fp256 ------------ 
    cout << " === implementation for Fp256 === " << endl;
    mpz_class arr[10001];
    for (int i = 1; i <= 10000; i++){arr[i] = r.get_z_range(p256);}    

    clock_t start, end; long double endtime;
    start = clock(); //程序开始计时
    for (int i = 1; i <= 10000; i++)
    {
        Hush_for_Huff_f(result, a256_1, b256_1, c256_1, d256_1, e256_1, arr[i], p256, exp256_f);
    }
    end = clock(); //程序结束用时
    endtime = (long double)(end - start) / CLOCKS_PER_SEC;
    cout << "    Total time:" << endtime * 1000 << "ms" << endl; //ms为单位

    return 0;
}

编译运行:

[hanss@Mint]$ g++ test_asm.cpp -o test -lgmp -lgmpxx
[hanss@Mint]$ ./test
 === implementation for Fp256 === 
    Total time:35.018ms

---

题外话:通用的计数程序写法(封装)

假设我们要对如下函数计数(可以发现它们的原型一致,即输入参数类型高度统一):

void Hush_for_Huff_f_s(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, mpz_class u, mpz_class p, mpz_class exp){...}
void Hush_for_Huff_f(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, mpz_class u, mpz_class p, mpz_class exp){...}
void constant_Hush_for_Huff_f_s(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, mpz_class u, mpz_class p, mpz_class exp, mpz_class exp_l, mpz_class exp_i){...}
void constant_Hush_for_Huff_f(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, mpz_class u, mpz_class p, mpz_class exp){...}

那么我们可以使用如下方式来封装:

// 封装
typedef void (*Hush_for_Huff)(mpz_class result[2], mpz_class a1, mpz_class b1, mpz_class c1, mpz_class d1, mpz_class e1, mpz_class u, mpz_class p, mpz_class exp, mpz_class exp_l, mpz_class exp_i);
...
// 调用
double timing(Hush_for_Huff FUNCTION, mpz_class result[2], mpz_class &a1, mpz_class &b1, mpz_class &c1, mpz_class &d1, mpz_class &e1, mpz_class u[], mpz_class &p, mpz_class &exp, mpz_class &exp_l, mpz_class &exp_i)
{
    ...
    FUNCTION(result, a1, b1, c1, d1, e1, u[i], p, exp, exp_l, exp_i);
    ...
}