В начало

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

Предлагается набор функций для работы с небольшими (до 2⁶⁴, для криптографии - недостаточно) простыми числами. Они организованы в виде статических членов класса 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;

Большая часть функций очевидна и не нуждается в комментариях. Поэтому прокомментируем лишь некоторые.

isDiv357_11_13:

  static bool isDiv357_11_13(UINT64 n); // n is div by 3,5,7,11,13?

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,...

Для совсем небольших простых чисел (до ~8000) все проверки проводятся по заранее заготовленным таблицам.

isPrime:

  static bool isPrime( UINT64 n); // is prime(n)?

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

Кроме того, метод не работает на квадратах простых чисел, поэтому их рассматриваем отдельно:

  UINT64 sn = (UINT64) sqrt ((double)n);
  if( sn*sn==n ) return sn;

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-битовых чисел: