В начало

Небольшие простые числа на C++

Предлагается набор функций для работы с небольшими (до 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




free counters