Предлагается набор функций для работы с небольшими (до 264,
для криптографии - недостаточно) простыми числами. Они организованы в виде статических
членов класса z64 (динамических членов у класса нет).
Модуль состоит из двух файлов:
z64.h - заголовочный файл
z64.cpp - реализация.
Основные методы класса:
typedef unsigned long long UINT64; typedef long long INT64; class z64 { ... public: static bool isFullSq (UINT64 n); // n is full square? static bool isDiv357_11_13(UINT64 n); // n is div by 3,5,7,11,13? static UINT64 gcd (UINT64 a, UINT64 b); // GCD(a,b) static UINT64 mul_mod (UINT64 a, UINT64 b, UINT64 с);// (a*b)%c static UINT64 pow_mod (UINT64 a, UINT64 b, UINT64 c);// (a^b)%c static bool MillerRabin(UINT64 a, UINT64 n); // MR primality test for n with base a static void euclid (UINT64 a, UINT64 b, INT64 &x, INT64 &y, UINT64 &d); // find x,y,d: a*x+b*y=d=gcd(a,b) static UINT64 inverse (UINT64 a, UINT64 n); // 1/a mod n static bool chineese (UINT64 a1, UINT64 b1, UINT64 a2, UINT64 b2, UINT64 &a3, UINT64 &b3); // cheneese theorem // sieve of Eratosthenes static void Efill(UINT64 mN); // fill sieve for primes till mN static void Efree(); // free memory static int Esave(char *fname); // return error code static int Eload(char *fname); // return error code // static UINT64 maxN(); // return maximal N in sieve static UINT64 nPrimes(); // total number of primes in the sieve (including 2) static UINT64 maxP(); // return maximal prime in the sieve // Primes static bool isPrime( UINT64 n); // is prime(n)? static UINT64 p( UINT64 i); // i-th prime, Maple: ithprime(i), p(0)=0, p(1)=2, p(2)=3,... static UINT64 Pi( UINT64 n); // the number of primes <= n, Pi(0)=Pi(1)=0, Pi(2)=1, Pi(3)=Pi(4)=2,... static UINT64 NextPrime( UINT64 n); // prime p>n, NextPrime(3)=5, NextPrieme(4)=5, NextPrime(0)=NextPrime(1)=2, NextPrime(2)=3 static UINT64 PrevPrime( UINT64 n); // prime p<n, PrevPrime(0)=PrevPrime(1)=PrevPrime(2)=0, PrevPrime(4)=3 // static UINT64 primeFactor( UINT64 n);// prime divisor or 1, if not found static UINT64 order( UINT64 a, UINT64 n);// order of a by modulo prime n static int jacobi(UINT64 a, UINT64 n);// Jacobi symbol static UINT64 msqrt( UINT64 a, UINT64 p);// sqrt(a) mod p, p must be prime! static char * factorString(UINT64 n); // string with factorisation: 180200 -> "2^3*5^2*17*53" static void toString(char *s, UINT64 n, int w=0); // s <- to String of width w static char *toString(UINT64 n, int w=0); // return to String of width w };Настройка пакета на различные процессоры и компиляторы сводится к корректировке двух строк в заголовочном файле:
typedef unsigned long long UINT64; typedef long long INT64;В ней надо указать название типа 64-битовых целых знаковых и беззнаковых чисел.
Большая часть функций очевидна и не нуждается в комментариях. Поэтому прокомментируем лишь некоторые.
isDiv357_11_13:
static bool isDiv357_11_13(UINT64 n); // n is div by 3,5,7,11,13?- проверка делимости на 3,5,7,11 и 13 за одно деление. Находим остаток от деления n на 3*5*7*11*13=15015 и затем по таблице (1877 байт) проверяем, будет ли остаток взаимно прост с 15015.
sieve of Eratosthenes: Все простые числа до 4 млрд. занимают 256 Мбайт памяти и вычисляются за 72 сек (AMD Athlon II X2 240, 2.8 GHz).
В решете хранятся только нечетные числа, по одному биту на каждое. Поэтому N байтов достаточно для храненения простых чисел до 16*N+1. То есть для хранения всех простых до 1 млн. требуется 63К байт. В одном гигабайте помещаются простые числа до 16 миллиардов.
С помощью решета (в пределах до maxP()) легко реализуются функции isPrime, NextPrime, PrevPrime. Однако для фукнций p() (ithprime) и Pi() требуются дополнительные данные. Одновременно с построением решета, запоминаются в отдельном массиве значения Pi(n) (количество простых, не больших n) для каждого n, кратного 1024. Таким образом, для хранения самого решета до числа N требуется N/16 байт, и ешё N/256 байт для этой дополнительной таблицы. С её помощью быстро работают функции
static UINT64 p( UINT64 i); // i-th prime, Maple: ithprime(i), p(0)=0, p(1)=2, p(2)=3,... static UINT64 Pi( UINT64 n); // the number of primes <= n, Pi(0)=Pi(1)=0, Pi(2)=1, Pi(3)=Pi(4)=2,...Эти функции взаимно-обратны: z64::Pi(z64::p(i)) == i.
Для совсем небольших простых чисел (до ~8000) все проверки проводятся по заранее заготовленным таблицам.
isPrime:
static bool isPrime( UINT64 n); // is prime(n)?В пределах имеющегося решета Эратосфена простота числа проверяется с его помощью. Если число n меньше 232 и не попадает в решето, то проводится проверка с помощью метода Миллера-Рабина по основаниям 2 и 3 и затем проверяются по таблице все 103 составных числа (unsigned BCompos[103]), меньшие 232, на которых ошибается метод. Таким образом, для чисел меньших 232 ответ всегда точный (а не вероятностный). Для больших чисел пока просто проводится 12 тестов Миллера-Рабина по основаниям 2,3,...,37.
mul_mod
Это - ключевая функция пакета. Реализация ее на ассемблере ускоряет работу
более, чем на порядок. Ассемблерная реализация проверена лишь для gcc-64 под Linux.
UINT64 z64::mul_mod(UINT64 a, UINT64 b, UINT64 c) // (a*b)%c { #if defined(_MSC_VER) || defined( _WIN32 ) UINT64 res=0; if(a>c) a %= c; while(b){ if(b&1) {res+=a; if(res>c) res-=c; b--; } else { a<<=1; if(a >c) a -=c; b>>=1;} } return res; #else UINT64 d; // to hold the result of a*b mod c // calculates a*b mod c, stores result in d asm ("mov %1, %%rax;" // put a into rax "mul %2;" // mul a*b -> rdx:rax "div %3;" // (a*b)/c -> quot in rax remainder in rdx "mov %%rdx, %0;" // store result in d :"=r"(d) // output :"r"(a), "r"(b), "r"(c) // input :"%rax", "%rdx" // clobbered registers ); return d; #endif } // z64::mul_mod
primeFactor:
static bool primeFactor( UINT64 n); // // prime divisor or 1, if not foundПо результатам некоторых исследований и большого числа экспериментов в наших пределах (до 264) из элементарных методов наиболее эффективен \rho-метод Полларда, все остальные не стоит и реализовывать. Если все простые множители у числа велики (несколько миллиардов), то методу может потребоваться до 2-3 сотен тысяч шагов. Но остальные методы - все равно ещё медленнее. Но, например, разложение числа p-1 на множители для простого p (это используется в функции order) находится без труда.
Кроме того, метод не работает на квадратах простых чисел, поэтому их рассматриваем отдельно:
UINT64 sn = (UINT64) sqrt ((double)n); if( sn*sn==n ) return sn;Эта проверка на числах, меньших 264 всегда является точной (проверено!).
factorString:
static char *factorString(UINT64 n); // string with factorisation: 180200 -> "2^3*5^2*17*53"Эта функция использует всегда одну и ту же статическую область памяти для формирования строки. То есть возвращаемый указатель - всегда один и тот же. Память не надо не выделять, ни освобождать. Надо только не забывать, что при повторном вызове этой функции, она снова использует ту же самую память. Если требуется преобразовать в строки сразу несколько чисел, то можно использовать функцию
static void toString(char *s, UINT64 n, int w=0); // s <- to String of width wНо о выделении достаточной памяти тогда уж придется заботится самому.
Пример использования:
#include "z64.h" ... z64::Efill(1<<24); // fill the sieve till 2^24 z64::Esave("16M.sieve"); // save the sieve UINT64 p = z64::p(1000000); // p = 15485863 UINT64 i = z64::Pi(p); // i = 1000000 UINT64 p2= z64::p( 700000); // p = 10570841 UINT64 n = p*p2; // n = 163 698 595 520 783 bool r = z64::isPrime(n); // r = false UINT64 f = z64::primeFactor(n); // f = 10570841
В основном, модуль разработан для учебных целей, и скорость работы не являлась главной целью. Однако при размере чисел в 64 бита, все равно вычисления получаются более быстрыми,чем при использовании библиотеки GMP. Например операция mul_mod на процессоре AMD Athlon II X2 240, 2.8 GHz занимает для 63-битовых чисел:
разрядность | язык | время, mks |
---|---|---|
32 | C++ | 75 |
64 | C++ | 31 |
64 | asm | 2.7 |
64 | GMP | 7 |