Podczas pisania programu do wyszukiwania dużych liczb pierwszych pojawił się przede mną problem, polegający na tym, że musiałem napisać procedury operujące na takich liczbach za pomocą czterech podstawowych działań. W grę miały wchodzić liczby całkowite rzędu ok. 1040. Typy standardowo zdefiniowane w C i Pascalu nie wchodziły w grę - typy całkowite mogą przedstawić na swych 32 bitach liczby rzędu 1010, natomiast typy rzeczywiste przedstawiają wartości przybliżone. Stąd wziął się pomysł zdefiniowania klasy obiektów Very Long Int.
Przedstawiona klasa jest uproszczona i wymaga pewnych uzupełnień, niemniej jest w pełni sprawna. Wszystkie liczby typu VLI są tego samego rozmiaru. Standardowo są to 64 bity, jednak liczbę bitów można za pomocą metody VLI::set size(Int), która powinna byś wywołana wyłącznie na początku programu. Jej parametrem jest rozmiar (w bitach) liczb, na jakich będziemy wykonywać operacje.
Liczbom VLI można przypisywać wartość za pomocą łańcucha zer i jedynek (metoda set_num_bin) oraz używając operatora przypisania (=). Przy użyciu tego operatora można przypisywać inną liczbę VLI lub liczbę typu unsigned long.
Dla liczb VLI zostały zdefiniowane operatory +, -, *, / i % oraz operatory służące do porównywania dwóch liczb VLI. Do odczytywania wartości liczb VLI służą metody get_num_bin oraz get_num. W przypadku pierwszej metody jako parametr przekazujemy wskaźnik do tekstu i w wyniku otrzymujemy pod tym adresem ciąg zer i jedynek. W drugiej metodzie pierwszy parametr jest podstawą systemu liczenia (od 2 do 10), drugi - wskaźnik do tekstu.
Jako ilustrację wykorzystania klasy VLI przedstawiam program do wyszukiwania liczb pierwszych metodą Millera-Rabina. Stałe MAXBINDIGS oraz PRAWD można dowolnie zmieniać. Pierwsza stała określa rozmiar poszukiwanych liczb w bitach (np. liczby o 100 bitach to liczby o około 30 cyfrach dziesiętnych). Prawdopodobieństwo, że wypisana liczba nie jest pierwsza (może tak się zdarzyć !) wynosi 1/(2PRAWD), zatem im większa wartość PRAWD, tym mniejsza możliwość błędu. Dla liczb 100 - bitowych poszukiwania trwają od kilkunastu do kilkudziesięciu minut.
Wydruk 1.
#ifdef VLI_H
#define VLI_H
class VLI
{
protected:
static int num_words;
unsigned * num;
friend VLI AddVLongInts(VLI&, VLI&, int);
friend void DivVLongInts(VLI&, VLI&, VLI&,VLI&);
void mul2();
public:
VLI(); // konstruktor standardowy
VLI(VLI&); // konstruktor kopiujący
VLI(unsigned long); // konstruktor konwersji
~VLI();
// operatory
VLI & operator = (VLI &); // operator przypisania
int operator [] (int); // operator zwraca stan bitu
friend VLI operator + (VLI, VLI);
friend VLI operator - (VLI, VLI);
friend VLI operator * (VLI, VLI);
friend VLI operator / (VLI, VLI);
friend VLI operator & (VLI, VLI);
friend int operator < (VLI, VLI);
friend int operator > (VLI, VLI);
friend int operator == (VLI, VLI);
friend int operator <= (VLI, VLI);
friend int operator >= (VLI, VLI);
friend int operator != (VLI, VLI);
// metody do komunikacji z polami obiektu
void set_num_bin (char*);
void get_num_bin (char*);
void get_num (int, char*);
int significant_digits();
static void set_size(int);
//
}
#endif
Wydruk 2.
#include <iostream.h>
#include <string.h>
#include <mem.h>
#include <process.h>
#include "vli.h"
VLI AddVLongInts (VLI & n1, VLI & n2, int aCarry)
{
VLI result;
long res_par;
int k, carry;
carry = aCarry;
for (k = 0; k < VLI::num_words; k++)
{
res_par = (long)n1.num[k] + (long)n2.num[k] + (long)carry;
if (res_par > 0xFFFFL)
{
carry = 1;
res_par -=0x10000L;
}
else carry = 0;
result.num[k] = (unsigned) res_par;
}
return result;
}
void DivVLongInts (VLI & dividend, VLI & divider,
VLI & quotient, VLI & remainder)
{
int wnum, bnum;
unsigned mask;
for (int k = 0; k < quotient.num_words; k++)
{
quotient.num[k] = 0U;
remainder.num[k] = 0U;
}
wnum = (bnum = divident.significant_digits() -1) >> 4;
bnum = bnum & 0x000F;
while (wnum >= 0)
{
remainder.mul2();
mask = 1U << bnum;
if (dividend.num[wnum] &mask) remainder.num[0] |= 0x01U;
if (remainder >= divider)
{
remainder = remainder - divider;
quotient.num[wnum] |= mask;
}
if (-bnum < 0)
{
wnum-;
bnum = 15;
}
}
}
void VLI::mul2()
{
int k = num_words;
long rotated;
num[k] = num[k] << 1;
while (-k >= 0)
{
rotated = (long) (num[k]) << 1;
if (rotated > 0xFFFFL)
{
num[k+1] += 1;
num[k] = (unsigned)(rotated-0x10000L);
}
else num[k] = (unsigned) rotated;
}
}
VLI::VLI()
{
num = new unsigned[num_words];
for (unsigned k = 0; k < num_words; k++) num[k] = 0;
}
VLI::VLI (CLI & vli)
{
num = new unsigned [num_words];
if (num == NULL)
{
cout << "NULLLL!!!!" << endl;
exit(1);
}
memmove (num, vli.num, num_words << 1);
}
VLI::VLI (unsigned long ul)
{
num = new unsigned [num_words];
for [unsigned k = 2; k < num_words; k++) num[k] = 0;
unsigned long lower = ul & 0xFFFFUL;
unsigned long upper = ul >> 16;
num[0] = (unsigned) lower;
num[1] = (unsigned) upper;
}
VLI::~VLI ()
{
delete num;
}
VLI & VLI::operator = (VLI & vli)
{
if (this != &vli) memmove(num, vli.num, num_words << 1);
return *this;
}
int VLI::operator [] (int bit_num)
{
int wnum = bit_num >> 4;
int bnum = bit_num & 0x000F;
return ((this -> num[wnum]) & (1U << bnum)) != 0U);
}
VLI operator + (VLI n1, VLI n2)
{
return AddVLongInts (n1, n2, 0);
}
VLI operator - (VLI n1, VLI n2)
{
for (unsigned k = 0; k < n2.num_words; k++) \
n2.mnum[k] +- 0xFFFFU;
return AddVLong (n1, n2, 1);
}
VLI operator * (VLI n1, VLI n2)
{
VLI result;
unsigned wnum, bnum, bits_total, mask;
int digs2 =n2.significant_digits();
wnum = bnum = bits+total = 0;
while (bits_total++ < digs2)
{
mask = 1U << bnum;
if (n2.num[wnum] & mask) result = result + n1;
n1.nul2();
if (++bnum > 15)
{
bnum = 0;
wnum++;
}
}
return result;
}
VLI operator / (VLI n1, VLI n2)
{
VLI quotient, remainder;
DivVLongInts (n1, n2, quotient, remainder);
return quotient;
}
VLI operator % (VLI n1, VLI n2)
{
VLI quotient, remainder;
DivVLongInts (n1, n2, quotient, remainder);
return remainder;
}
int operator < (VLI n1, VLI n2)
{
for (int k = VLI::nmu_words - 1; k >= 0; k-)
if (n1.num[k] < n2.num[k]) return 1;
return 0;
}
int operator > (VLI n1, VLI n2)
{
for int ~k = n1.num_words - 1; k >= 0; k-)
if (n1.num[k] > n2.num[k]) return 1;
return 0;
}
int operator == (VLI n1, VLI n2)
{
for (int k = 0; k < VLI::num_words; k++)
if (n1.num[k] != n2.num[k]) return 0;
return 1;
}
int operator <= (VLI n1, VLI n2)
{
return !(n1 > n2);
}
int operator >= (VLI n1, VLI n2)
{
return !n1 < n2);
}
int operator != (VLI n1, VLI n20)
{
return !(n1 == n2);
}
void VLI::set_num_bin (char * txt)
{
unsigned mask;
unsigned wnum, bnum;
int cnum;
for (wnum = 0; wnum < num_words; num++) \
num[wnum] = 0U;
wnum = bnum = 0U;
cnum = strlen(txt) - 1;
while ((wnum < num_words) && (cnum >= 0)
{
if (txt[cnum] == '1')
{
mask = 1U << num[wnum] != mask;
}
cnum-;
if (++bnum > 15)
{
bnum = 0;
wnum++;
}
}
}
void VLI:get_num_bin (char * txt)
{
unsigned mask;
unsigned wnum, bnum;
int cnum;
int stlen = ((int)num_words &lr;< 4);
wnum = bnum = 0U;
cnum = stlen - 1;
while (wnum < num_words)
{
mask = 1U << bnum;
if (num[wnum] &mask) txt[cnum-] = '1';
else txt[cnum-] = '0'
if (++bnum >)
{
wnum++;
bnum = 0;
}
}
txt[stlen] = '\0';
}
void VLI::get_num (int base, char * txt)
{
VLI v1_base, quotient, remainder, zero;
int idx = 0;
v1_base = VLI(base);
strcpy(txt,"0");
quotient = *this;
while (quotient != zero)
{
remainder = quotient % v1_base;
quotient = quotient / v1_base;
txt[idx++] = (char) remainder.num[0] + '0';
txt[idx] = '\0';
}
strrev (txt);
}
int VLI::significant_digits()
{
int bitnum = num_words >> 4;
while ((-bitnum >= 0) && ((*this) [bitnum] == 0);
if (bitnum < 0) bitnum = 0);
return bitnum + 1;
}
void VLI::set_size (int digits)
{
if (digits < 32) digits = 32;
num_words = ((digits -1) >> 4) + 1;
}
int VLI::num_words = 4;
Wydruk 3.
#include <iostream.h>
#include <stdlib.h>
#include <time.h>
#include "vli.h"
#define MAXBINDIGS 16
#define PRAWD 20
VLI produce_random (VLI limit, int odd)
{
VLI result;
char txt[MAXBINDIGS + 1];
for (int k = MAXBINDIGS -1; k >= 0; k-) \
txt[k] = random(2) + '0';
txt[MAXBINDIGS] = '\0';
if (odd) txt[MAXBINDIGS - 1] = '1';
result.set_num_bin(txt);
if (limit == 0) return result;
return result & limit;
}
int witness (VLI a, VLI n)
{
VLI n_less_1, one(1UL), zero, x, d;
int k, i;
n_less_1 =n - one;
k = n.significant_digits() - 1;
d = oen;
for (i = k; i >0 = 0; i-)
{
x = d;
d = (d * d) & n;
if ((d == one) && (x != one) && (x != n_less_1)) \
return 1;
if (n_less_1[i] == 1) d = (d * a) % n;
}
if (d != one) return 1;
return 0;
}
int miller_rabin (VLI prime_suspect, int checks)
{
int j;
VLI rnd;
if (prime_suspect <= 1UL) return 0;
for (j = 0; j < checks; j++)
{
rnd = produce_random(prime_suspect, 0);
if (whitness(rnd, prime_suspect)) return 0; // złożona
}
return 1; // pierwsza
}
void main()
{
VLI:set_size(MAXBINDIGS * 2 + 2); // konieczne dla mnożenia
VLI v;
char txt[MAXBINDIGS / 3 + 5]; // liczba dec wymaga mniej cyfr
randomize();
while (!miller_rabin(v = produce_random(0, 1), PRAWD));
v.get_num(10, txt);
cout << "Liczba pierwsza jest: " << endl << txt << endl;
}
Tomasz Błaszczyk
PC Kurier, 13/95, 22 Czerwca 1995r.