メルセンヌ素数2^{521}-1を法とする場合の乗算の実行時間比較

mulMod521.cpp

#include <iostream>
#include <time.h>

#include <gmpxx.h>

#define SIZE 9

void mod2(mp_limb_t *z, const mp_limb_t *XY, const mp_limb_t *p, mp_limb_t *q) {
    mpn_tdiv_qr(q, z, 0, (const mp_limb_t *)XY, SIZE*2, (const mp_limb_t*)p, SIZE);
}

void mod3(mp_limb_t *z, const mp_limb_t *XY, const mp_limb_t *p, mp_limb_t *t, mp_limb_t *s) {
    // (T + (T mod R)*N) / R
    mpn_zero(t, SIZE*2);
    mpn_zero(s, SIZE*2);
    mpn_and_n(t, XY, p, SIZE); // T mod R
    mpn_mul_n(s, (const mp_limb_t*)t, p, SIZE);
    mpn_add_n(t, (const mp_limb_t*)s, XY, SIZE*2); // (T + (T mod R)*N)

    mpn_rshift(t, (const mp_limb_t*)t, SIZE*2, 9);
    for (size_t i = 0; i < SIZE; i++) { // (T + (T mod R)*N) / R
        t[i] = t[i+8];
    }
    if (mpn_cmp((const mp_limb_t*)t, p, SIZE) >= 0) {
        mpn_sub_n(t, (const mp_limb_t*)t, p, SIZE);
    }
    mpn_copyi(z, (const mp_limb_t*)t, SIZE);
}

void mod(mp_limb_t *z, const mp_limb_t *XY, const mp_limb_t *p, mp_limb_t *t, mp_limb_t *s) {
    // (T + (T mod R)*N) / R
    mpn_zero(t, SIZE*2);
    mpn_zero(s, SIZE*2);
    mpn_and_n(t, XY, p, SIZE); // T mod R
    for (size_t i = 0; i < SIZE; i++) {
        s[i+8] = t[i];
    }
    mpn_lshift(s, (const mp_limb_t*)s, SIZE*2, 9);
    mpn_sub_n(s, (const mp_limb_t*)s, (const mp_limb_t*)t, SIZE*2); // (T mod R)*(N+1) - (T mod R)

    mpn_add_n(t, (const mp_limb_t*)s, XY, SIZE*2); // (T + (T mod R)*N)

    mpn_rshift(t, (const mp_limb_t*)t, SIZE*2, 9);
    for (size_t i = 0; i < SIZE; i++) { // (T + (T mod R)*N) / R
        t[i] = t[i+8];
    }
    if (mpn_cmp((const mp_limb_t*)t, p, SIZE) >= 0) {
        mpn_sub_n(t, (const mp_limb_t*)t, p, SIZE);
    }
    mpn_copyi(z, (const mp_limb_t*)t, SIZE);
}

void mulMod(mpz_class &z, const mpz_class &x, const mpz_class &y, const mpz_class &p) {
    z = x*y % p;
}

int main() {
    mp_limb_t m[SIZE] = {
        0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff,
        0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff,
        0x1ff,
    };
    mp_limb_t a[SIZE] = {
        0x2137201477420138L, 0x6313232130472104L, 0x3639164821638291L, 0x8371372917431694L,
        0x1739278732747374L, 0x1732913789724789L, 0x7296392641693692L, 0x32dead13L,
    };
    mp_limb_t b[SIZE] = {
        0x8032810382412080L, 0x2427104721732301L, 0x3281037402174017L, 0x7201839284217487L,
        0x1737201372138214L, 0x3747217427103782L, 0x13729dee74L,
    };

    mp_limb_t ab[SIZE*2] = {0};
    mp_limb_t c[SIZE];

    mp_limb_t tp [SIZE*2] = {0};
    mp_limb_t sp [SIZE*2] = {0};

    const int N = 1000000;
    time_t begin, end;

    begin = clock();
    for (int i = 0; i < N; i++) {
        mpn_mul_n(ab, (const mp_limb_t*)a, (const mp_limb_t*)b, SIZE);
        mod2(c, (const mp_limb_t*)ab, (const mp_limb_t*)m, sp);
    }
    end = clock();
    printf("[+] mulMod tdiv mpn: \t%f usec\n", double(end-begin)*1e6/CLOCKS_PER_SEC/N);

    begin = clock();
    for (int i = 0; i < N; i++) {
        mpn_mul_n(ab, (const mp_limb_t*)a, (const mp_limb_t*)b, SIZE);
        mod(c, (const mp_limb_t*)ab, (const mp_limb_t*)m, tp, sp);
    }
    end = clock();
    printf("[+] mulMod MR1 mpn: \t%f usec\n", double(end-begin)*1e6/CLOCKS_PER_SEC/N);

    begin = clock();
    for (int i = 0; i < N; i++) {
        mpn_mul_n(ab, (const mp_limb_t*)a, (const mp_limb_t*)b, SIZE);
        mod3(c, (const mp_limb_t*)ab, (const mp_limb_t*)m, tp, sp);
    }
    end = clock();
    printf("[+] mulMod MR2 mpn: \t%f usec\n", double(end-begin)*1e6/CLOCKS_PER_SEC/N);


    mpz_class p("1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 16);
    mpz_class x("32DEAD137296392641693692173291378972478917392787327473748371372917431694363916482163829163132321304721042137201477420138", 16);
    mpz_class y("13729DEE74374721742710378217372013721382147201839284217487328103740217401724271047217323018032810382412080", 16);
    mpz_class z;
    begin = clock();
    for (int i = 0; i < N; i++) {
        mulMod(z, x, y, p);
    }
    end = clock();
    printf("[+] mulMod mpz: \t%f usec\n", double(end-begin)*1e6/CLOCKS_PER_SEC/N);
}

mulMod_log

$ cat /proc/cpuinfo | grep "model name"
model name	: Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz


$ g++ mul521.cpp  -lgmpxx -lgmp -Ofast && ./a.out

"""
[+] mulMod tdiv mpn: 	0.187029 usec
[+] mulMod MR1 mpn: 	0.098798 usec
[+] mulMod MR2 mpn: 	0.129871 usec
[+] mulMod mpz: 	0.178452 usec
"""