Для алгоритма шифрования Эль-Гамаля надо найти большие простые числа порядка 2^24-2^32

Написал простую прогу.
Основана на малой теореме Ферма

до числа 80476 она работает отлично, однако, при базе =5, число 80476 выдается как простое. дальше количество простых чисел на одном и том же интервале зависит от значения базы.

Если у Вас есть алгоритмы нахождения таких простых чисел или Вы видите ошибку моей программе напишите пожалуйста.

Мих, Falcon17@pisem.net

Заранее благодарен.

текст
/*******/
#include <stdio.h>
#include <conio.h>

/*

основана на тмалой теореме ферма:
если х**(р-1) % р == 1, то р — простое х=3,5,7
Здесь 2 функции:
powmod
возвр r = (a**b)%p
multrest
return (a*b) %p

*/

unsigned long mulrest(unsigned long a,unsigned long b,unsigned long c)
{
unsigned long A[2], B[2], X[3], x0r, x1r;
A[0] = a>>16; //разбиваем числа a b надвое
A[1] = a — A[0]*65536;//32768;
B[0] = b>>16;
B[1] = b — B[0]*65536;//32768;
X[0] = A[0]*B[0];
X[1] = A[1]*B[0] + A[0]*B[1];
X[2] = A[1]*B[1];

x0r = (4294967295 % c + 1)%c;
x1r = 65536 % c;

X[0] = (X[0]%c) * x0r;
X[1] = (X[1]%c) * x1r;
X[2] = X[2]%c;

return (X[0]+X[1]+X[2])%c;
}

unsigned long int powmod(unsigned long int a,unsigned long int b,unsigned long p)
{
unsigned long int r=1, t[32];
int i=0;
//быстрое возведение в степень
t[0]=a; // заполним массив t числами a, a*a %p, ...
for (i=1;i<32;i++)
t[i]=mulrest(t[i-1],t[i-1],p);

for (i=1;i<32;i++)
if (((b>>(i-1))%2)) // представим b в двоичном коде, где 1
r=mulrest(r,t[i-1],p); //r умножаем на t[i-1]

return r;
}

void main()
{
clrscr();
unsigned long int i=0,q=0,j;
unsigned long int d,e,s,m,base = 3;
d=1; //от какого числа
e=100; //диапазон
for(i=0;i<e;i++)
{
m=base;
s=d++;
m=powmod(m,s-1,s);
if(m==1) {printf("%7lu|%9lu ",++q,s);
if((q+1)%2)printf("\n");}
}
getch();
}