Год назад я писал о простых числах Софи Жермен.

Вкратце напомню: p — простое число Софи Жермен, если q = 2p+1, тоже простое число. Простые числа Софи Жермен применяются в криптографии (в частности, в протоколе обмена ключами Диффи–Хеллмана–Меркле).

Ранее я использовал вероятностный метод Рабина–Миллера для проверки обоих чисел (p и q) на простоту. Вероятностная проверка (вкупе с другими тестами, например, решетом Эратосфена и теста на псевдопростоту, основанного на малой теореме Ферма) — не самый быстрый (и не всегда точный) метод.

Как оказалось, есть более простая и эффективная проверка, которую открыл Генри Лифшиц (Henry Lifchitz).

Суть проверки сводится к тому, что если p≥5 является простым числом, то q=2p+1 является простым числом тогда и только тогда, когда q является делителем числа 3p-1.

Среди преимуществ такого подхода метода можно отметить его простоту. Среди недостатков — при реализации алгоритма с помощью GMP сразу сталкиваешься с тем, что в библиотеке нет функции mpz_pow(res, x, y), такой, чтобы аргументы x и y имели тип mpz_t. Судя по всему, и не будет. В принципе, логично: возведение тройки в степень числа с 4096 разрядами съело больше 4 ГиБ памяти и я снял задачу прежде, чем OOM Killer убил какой-нибудь процесс.

Кому интересно поэкспериментировать, привожу программу, находящую простые числа Софи Жермен до 1,000,000:

[-]
Download sophie1.c
#include <stdlib.h>
#include <gmp.h>

static void sj_mpz_pow(mpz_ptr y, mpz_srcptr x, mpz_srcptr n)
{
    mpz_t i;
    mpz_t z;

    mpz_set_ui(y, 1);
    mpz_set(i, n);
    mpz_set(z, x);

    do {
        if (0 != mpz_odd_p(i)) {
            mpz_mul(y, y, z);
        }

        mpz_mul(z, z, z);
        mpz_divexact_ui(i, i, 2);
    } while (0 != mpz_cmp_ui(i, 0));
}

int main(void)
{
    mpz_t x;
    mpz_t q;
    mpz_t r;
    mpz_t three;

    mpz_t three_powered;
    mpz_t last_power;
    mpz_t pow_diff;

    mpz_init_set_ui(three_powered, 3);
    mpz_init_set_ui(last_power, 1);
    mpz_init(pow_diff);

    mpz_init_set_ui(x, 5);
    mpz_init_set_ui(three, 3);
    mpz_init(r);
    mpz_init(q);

    do {
        mpz_add(q, x, x);
        mpz_add_ui(q, q, 1); /* q = 2x + 1 */

        mpz_sub(pow_diff, x, last_power);

        if (0 != mpz_fits_uint_p(pow_diff)) {
            unsigned int diff = mpz_get_ui(pow_diff);
            mpz_pow_ui(r, three, diff);
            mpz_mul(three_powered, three_powered, r);
        }
        else {
            sj_mpz_pow(r, three, pow_diff);
            mpz_mul(three_powered, three_powered, r);
        }

        mpz_set(last_power, x);
        mpz_sub_ui(r, three_powered, 1);

        if (mpz_divisible_p(r, q)) {
            gmp_printf("%Zd\n", x);
        }

        mpz_nextprime(x, x);
    } while (mpz_cmp_ui(x, 1000000) < 0);

    return EXIT_SUCCESS;
}

Для вычисления степени приведён (но реально не используется) алгоритм addition-chain exponentiation.

Это было решение в лоб, и оно не очень (или очень не) красиво. Если же вспомнить математику, то всё становится гораздо проще. Математику уже затем учить следует, что она ум в порядок приводит (М.В. Ломоносов).

Ведь на самом деле,

(a + b) mod c = a mod c + b mod c

Из этого следует, что

(3p-1) mod (2p+1) = 3p mod (2p+1) - 1 mod (2p+1)

Иными словами, если выполняется равенство

3p mod (2p+1) = 1,

то p — простое число Софи Жермен.

Что это даёт: мы можем использовать оптимизированный алгоритм для вычисления ab mod c.

В результате получим очень простую программу:

[-]
Download sophie2.c
#include <stdlib.h>
#include <gmp.h>

int main(void)
{
    mpz_t x;
    mpz_t q;
    mpz_t r;
    mpz_t three;

    mpz_init_set_ui(x, 5);
    mpz_init_set_ui(three, 3);
    mpz_init(r);
    mpz_init(q);

    do {
        mpz_add(q, x, x);
        mpz_add_ui(q, q, 1); /* q = 2x + 1 */

        mpz_powm(r, three, x, q);

        if (0 == mpz_cmp_ui(r, 1)) {
            gmp_printf("%Zd\n", x);
        }

        mpz_nextprime(x, x);
    } while (mpz_cmp_ui(x, 1000000) < 0);

    return EXIT_SUCCESS;
}

Такие вот пироги
А выигрыш в производительности очень заметен на действительно больших числах.

Добавить в закладки

Связанные записи

22
Апр
2009

Комментарии к статье «И снова о простых числах Софи Жермен»  »

К статье «И снова о простых числах Софи Жермен» комментариев пока нет. Не хотите ли стать первым?

Подписаться на RSS-ленту комментариев к статье «И снова о простых числах Софи Жермен» Trackback URL: http://blog.sjinks.org.ua/c-cpp/548-on-sophie-germain-primes-one-again/trackback/

Оставить комментарий к записи «И снова о простых числах Софи Жермен»

Вы можете использовать данные тэги: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

Оставляя комментарий, Вы выражаете своё согласие с Правилами комментирования.

Подписаться, не комментируя