Here is a C implementation for word-sized integers:
100:#include <stdio.h>
101:#include <fcntl.h>
102:#include <stdlib.h>
103:#include <stdint.h>
104:#include <inttypes.h>
105:#include <string.h>
106:#include <errno.h>
107:#include <unistd.h>
108:#include <sys/stat.h>
109:#include <stdarg.h>
110:
111:
112:void SYSe(const char *fmt,...)
113:{
114: va_list ap;
115: int eno;
116: eno= errno;
117: fprintf(stderr,"SYS error: ");
118: va_start(ap,fmt);
119: vfprintf(stderr,fmt,ap);
120: va_end(ap);
121: fprintf(stderr," %s\n", strerror(eno));
122: exit(1);
123:}
124:
125:void read_urandom(unsigned char *d, int len)
126:{
127: int fd,R;
128: fd= open("/dev/urandom", O_RDONLY);
129: if (fd<0) SYSe("open /dev/urandom");
130: R= read(fd, d, len);
131: if (R!=len) SYSe("read %d returned %d", len, R);
132: close(fd);
133:}
134:
135:uint32_t rand32(uint32_t lo, uint32_t hi)
136:{
137: uint32_t V;
138: while(1)
139: {
140: read_urandom((unsigned char*) &V, sizeof(V));
141: if (V>=lo && V<=hi) return V;
142: }
143:}
144:
145:uint32_t pow32n(uint32_t B, uint32_t P, uint32_t N)
146:{
147: uint64_t A,X;
148: A= 1;
149: X= B;
150: while(P)
151: {
152: if (P&1) A= A*X%N;
153: X= X*X%N;
154: P>>= 1;
155: }
156: return A;
157:}
158:
159:void separate_candidate(uint32_t *Rd, uint32_t *Rr, uint32_t N)
160:{
161: uint32_t r;
162: N--; N>>= 1; r= 1;
163: while((N!=0) && (N&1==0))
164: {
165: r++;
166: N>>= 1;
167: }
168: *Rd= N;
169: *Rr= r;
170:}
171:
172:int witness(int d, uint32_t r, uint32_t N)
173:{
174: uint64_t x;
175: int i;
176: x= pow32n(rand32(2,N-2), d, N);
177: if (x==1 || x==N-1) return 0;
178: for(i=0;i<r-1;i++)
179: {
180: x= x*x%N;
181: if (x==N-1) return 0;
182: }
183: return 1;
184:}
185:
186:int isprime(uint32_t N, int k)
187:{
188: uint32_t r,d;
189: int i;
190: separate_candidate(&d, &r, N);
191: for(i=0;i<k;i++)
192: if (witness(d,r,N)) return 0;
193: return 1;
194:}
195:
196:uint32_t genprime(int k)
197:{
198: uint32_t N;
199: while(1)
200: {
201: N= rand32(3,-1);
202: N |= 1;
203: N |= 1<<31; // let it be big enough..
204: if (isprime(N,k)) return N;
205: }
206:}
207:
208:
209:int main(int argc,char **argv)
210:{
211: printf("%" PRIu32 "\n", genprime(atoi(argv[1])));
212: return 0;
213:}
Let's discuss this code and improve it before making a multi-precision
version of it. First of all, let's look at how random numbers are
generated. I chose to ignore any random number that is not in the
given range because doing arithmetic manipulation to make it fit
may reduce the entropy of these numbers. There are only two calls
to the random number generator and the set of excluded numbers isn't
that big.
In (176), the excluded numbers are: 0,1,N-1 thru 2n-1 where 2n-1<N≤2n-1. In (201), we only have 0,1 and 2. For the latter, the set of excluded numbers is so small, it's practically impossible for a random number generator to return a sequence of such small numbers for an extended period of time.
For the first case, approximately half the numbers we get will be outside of the required range. Therefore, the loop in rand32() will run for 1.5 steps on average. The probability that the loop continues for t steps is 2-t. Theoretically, it could go on forever, but I'm willing to take that exceedingly small risk of a broken random number generator. Therefore, I will keep the random number generator as is.
Now let's look at the pow32n and witness functions. In these functions, d is not subjected to any arithmetic operation. It's only used for computing ad, where its bits are examined one by one. The same bits are already present in N, beyond the binary digit #r. There is no need to generate d at all. We can simply read the bits from N, but by skipping the first r bits.
500:uint32_t pow32nskip(uint32_t B, uint32_t N, uint32_t skip)
501:{
502: uint64_t A,X;
503: uint32_t bit;
504: A= 1;
505: X= B;
506: for(bit=skip;bit<32;bit++)
507: {
508: if (N&(1<<bit)) A= A*X%N;
509: X= X*X%N;
510: }
511: return A;
512:}
513:
514:uint32_t count_zeros(uint32_t N)
515:{
516: uint32_t r;
517: for(r=1;r<32;r++) if (N&(1<<r)) return r;
518: return 0; // impossible.
519:}
520:
521:int witness(uint32_t r, uint32_t N)
522:{
523: uint64_t x;
524: int i;
525: x= pow32nskip(rand32(2,N-2), N, r);
526: if (x==1 || x==N-1) return 0;
527: for(i=0;i<r-1;i++)
528: {
529: x= x*x%N;
530: if (x==N-1) return 0;
531: }
532: return 1;
533:}
534:
535:int isprime(uint32_t N, int k)
536:{
537: uint32_t r;
538: int i;
539: r= count_zeros(N);
540: for(i=0;i<k;i++)
541: if (witness(r,N)) return 0;
542: return 1;
543:}
We have a couple of improvements.
The only function I don't have in there is the powNskip function, which skips a given number of least-significant digits from the power.
I will need the constants 1 and N-1 in both Montgomery and regular form. The former for comparison in witness(526 and 530) and the latter for generating the candidate bases(525).
The number of trials for ensuring that a given number is prime with a given probability (k in the code) can be computed by using Appendix F.1 in [4]. The openssl source code contains a table of these values. I copied the numbers from there, I hope that there are no license issues with that.