К оглавлению

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

С.И.Хашин

Простых чисел, меньших 10¹¹ (100 млрд.) меньше, чем 2³², поэтому для хранения номеров таких простых чисел достаточно 4-х байт. На самом деле,

        π( 104.4 ⋅ 109 )= 4 291 624 115 = 232 - 3343181.

Заметим, что среди вычетов по модулю 30 ровно 8 взаимно просты с 30, это ±1, ±7, ±11, ±13. Поэтому удобно будет хранить в байте номер k признаки простоты чисел

        30k+1, 30k+7, 30k+11, 30k+13, 30k+17, 30k+19, 30k+23, 30k+29

pr[k]=номер простого числа m, если оно простое
      и -( наименьший простой делитель m) иначе,
      где m=2k+1.

i:     0   1   2   3   4   5   6   7   8   9  10  11  12  13  14 15  16  17  18  19
2*i+1: 1   3   5   7   9  11  13  15  17  19  21  23  25  27  29 31  33  35  37  39
pr[i]: 0   2   3   4  -3   5   6  -3   7   8  -3   9  -5  -3  10 11  -3  -5  12  -3
pi[i]: 0   2   3   5   7  11  13  17  19  23  29  31  37  41  43 47  53  59  61  67

Благодаря этому, для всех нечётных чисел, меньших R1 через одно обращение к таблице находится наибольший простой делитель числа.

Таким образом, в максимальном варианте, для хранения всех простых чисел, меньших 100 млрд. нам потребуется 3Гбайта оперативной памяти и ещё 400М байт для хранения наибольших простых делителей чисел, меньших 200 млн..
При современных объемах оперативной памяти — вполне умеренные требования!

Для реализации этого предлагается набор функций для работы с небольшими (до 2⁶⁴, для криптографии - недостаточно) простыми числами. Они организованы в виде статических членов класса z64 (динамических членов у класса нет). Модуль состоит из двух файлов:

z64.h - заголовочный файл
z64.cpp - реализация.

В модуле не используются никакие внешние библиотеки.
Единственная сложность – использование ассемблера в функции mul_mod.
Если вы не хотете с этим связываться, ассемблер можно отключить, все вычисления будут идти на С++. При этом функция mul_mod и связанные с ней (типа pow_mod) будут работать в 3-4 раза медленнее.

Основные методы класса:

typedef unsigned long long u64;
typedef          long long i64;

class z64 {
...
public:
    static u64 maxN();              // return maximal N
    static u64 maxP();              // return maximal prime in the sieve

    // Primes
    static bool isPrime( u64 n);    // is prime(n)?
    static u64 p ( u64 i);          // i-th prime,  Maple: ithprime(i), p(0)=0, p(1)=2, p(2)=3,...
    static u64 ip( u64 n);          // index of prime number, 0 if n is composite
    static u64 Pi( u64 n);          // the number of primes ≤ n, Pi(0)=Pi(1)=0, Pi(2)=1, Pi(3)=Pi(4)=2,...
    static u64 NextPrime( u64 n);   // prime p>n, NextPrime(0)=NextPrime(1)=2, NextPrime(2)=3, NextPrime(3)=NextPrime(4)=5.
    static u64 PrevPrime( u64 n);   // prime p<n, PrevPrime(0)=PrevPrime(1)=PrevPrime(2)=0, PrevPrime(4)=3

    // Factoring
    static u64  primeFactor( u64 n);// prime divisor or 1, if not found


    static bool    isFullSq (u64 n);               // n is full square?
    static u64     gcd      (u64 a, u64 b);        // GCD(a,b)
    static u64     mul_mod  (u64 a, u64 b, u64 n); // (a*b)%c
    static u64     pow_mod  (u64 a, u64 b, u64 c); // (a^b)%c
    static bool    MillerRabin(u64 a,u64 n);       // MR primality test for n with base a
    static void    euclid   (u64 a, u64 b, i64 &x, i64 &y, u64 &d); // find x,y,d: a*x+b*y=d=gcd(a,b)
    static u64     inverse  (u64 a, u64 n);        // 1/a mod n
    static bool    chineese (u64 mod1, u64 b1, u64 mod2, u64 b2, u64 &mod3, u64 &b3); // cheneese theorem

    // Rest
    static u64  order( u64 a, u64 n);   // order of a by modulo prime n
    static int  jacobi(u64 a, u64 n);   // Jacobi symbol
    static char *factorString(u64 n);   // string with factorisation: 180200 -> 2^3*5^2*17*53
    static u64  msqrt( u64 a, u64 p);   // sqrt(a) mod p, p must be prime!

    // Frobenius
    static void   Frobenius( u64 a, u64 b, i64 c, u64 k, u64 p, u64 &a1, u64 &b1);
                  // a1+b1*sqrt(c) = (a+b*sqrt(c))^k mod p
    static u64    FrobOrder( u64 a, u64 b, i64 c, u64 p);   // order(a+b*sqrt(c)) mod p
    static bool   FrobTest ( u64 a, u64 b, i64 c, u64 n);   // n is prime?
    static int    minC(u64 n);                              // minimal c: J(c/n)=-1;

    static void    toString(char *s, u64 n, int w=0);       // s <- to String of width w
    static char   *toString(u64 n, int w=0);                // return to String of width w

};

typedef unsigned long long u64;
typedef          long long i64;

mul_mod
Это - ключевая функция пакета. Реализация ее на ассемблере ускоряет работу более, чем на порядок. Ассемблерная реализация проверена лишь для gcc-64 под Linux.

u64  z64::mul_mod(u64 a, u64 b, u64 c)        // (a*b)%c
{
#if defined(_MSC_VER) || defined( _WIN32 )
    u64 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
  u64 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

Под Visual C++/64 функция также реализована на ассемблере в модуле mul_mod.asm:

; extern "C" u64     mul_mod_ams(u64 a, u64 b, u64 n);
; Project -> Properties -> Build Dependences -> Build Customizations.. -> пометить masm
; mul_mod.asm->Properties->(Platform=x64) Microsoft Macro Assembler->Preprocessor Definition->_WIN64
; команда mul умножает регистр rax на другой регистр
IFDEF _WIN64            ; В 32-разрядном режиме не исользуется
.CODE
mul_mod_asm PROC
;  rcx = a
;  rdx = b
;  r8  = n
    mov rax, rdx        ; put b into rax
    mul rcx             ; mul a*b -> rdx:rax
    div r8              ; (a*b)/n -> quot in rax remainder in rdx
    mov rax, rdx        ; Копирование RDX в RAX 
    ret                 ; Возврат в главную программу 
mul_mod_asm ENDP
ENDIF
END

primeFactor:

    static u64  primeFactor( u64 n);    // prime divisor or 1, if not found

По результатам некоторых исследований и большого числа экспериментов в наших пределах (до 2⁶⁴) из элементарных методов наиболее эффективен оказался ρ-метод Полларда, все остальные не стоит и реализовывать. Если все простые множители у числа велики (несколько миллиардов), то методу может потребоваться до 2-3 сотен тысяч шагов. Но остальные методы - все равно ещё медленнее. Но, например, разложение числа p-1 на множители для простого p (это используется в функции order) находится без труда.

factorString:

    static void    toString(char *s, u64 n, int w=0);       // s >- to String of width w
    static char   *toString(u64 n, int w=0);                // return to String of width w

Эта функция использует всегда одну и ту же статическую область памяти для формирования строки. То есть возвращаемый указатель - всегда один и тот же. Память не надо не выделять, ни освобождать. Надо только не забывать, что при повторном вызове этой функции, она снова использует ту же самую память. Если требуется преобразовать в строки сразу несколько чисел, то можно использовать функцию

  static void toString(char *s, u64 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
    u64 p = z64::p(1000000);     // p = 15485863
    u64 i = z64::Pi(p);          // i =  1000000
    u64 p2= z64::p( 700000);     // p = 10570841
    u64 n = p*p2;                // n = 163 698 595 520 783
    bool   r = z64::isPrime(n);  // r = false
    u64 f = z64::primeFactor(n); // f = 10570841

В основном, модуль разработан для учебных целей, и скорость работы не являлась главной целью. Однако при размере чисел в 64 бита, все равно вычисления получаются более быстрыми,чем при использовании библиотеки GMP. Например операция mul_mod на процессоре AMD Athlon II X2 240, 2.8 GHz занимает для 63-битовых чисел:

К оглавлению