1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "primes-lcg.h"
4 #include "primelist-lcg.h"
5
6 #define YES 1
7 #define NO 0
8 #define NPRIMES 1000
9
10 int primes[NPRIMES];
11
12 #ifdef __STDC__
init_prime(void)13 int init_prime(void)
14 #else
15 int init_prime()
16 #endif
17 {
18 int i, j, obtained = 0, isprime;
19
20 for(i=3; i < MINPRIME; i += 2)
21 {
22 isprime = YES;
23
24 for(j=0; j < obtained; j++)
25 if(i%primes[j] == 0)
26 {
27 isprime = NO;
28 break;
29 }
30 else if(primes[j]*primes[j] > i)
31 break;
32
33 if(isprime == YES)
34 {
35 primes[obtained] = i;
36 obtained++;
37 }
38 }
39
40 return obtained;
41 }
42
43
44
45
46 #ifdef __STDC__
getprime(int need,int * prime_array,int offset)47 int getprime(int need, int *prime_array, int offset)
48 #else
49 int getprime(need, prime_array, offset)
50 int need, *prime_array,offset;
51 #endif
52 {
53 static int initiallized = NO, num_prime;
54 int largest;
55 int i, isprime, index, obtained = 0;
56
57 if(need <= 0)
58 {
59 fprintf(stderr,"WARNING: Number of primes needed = %d < 1; None returned\n"
60 , need);
61 return 0;
62 }
63
64 if(offset < 0)
65 {
66 fprintf(stderr,"WARNING: Offset of prime = %d < 1; None returned\n"
67 , offset);
68 return 0;
69 }
70
71
72 if(offset+need-1<PRIMELISTSIZE1)
73 {
74 memcpy(prime_array,prime_list+offset,need*sizeof(int));
75 return need;
76 }
77
78 if(!initiallized)
79 {
80 num_prime = init_prime();
81
82
83 largest = MAXPRIME;
84 initiallized = YES;
85 }
86
87 if(offset > MAXPRIMEOFFSET)
88 {
89 fprintf(stderr,"WARNING: generator has branched maximum number of times;\nindependence of generators no longer guaranteed");
90 offset = offset % MAXPRIMEOFFSET;
91 }
92
93 if(offset < PRIMELISTSIZE1) /* search table for previous prime */
94 {
95 largest = prime_list[offset] + 2;
96 offset = 0;
97 }
98 else
99 {
100 index = (int) ((offset-PRIMELISTSIZE1+1)/STEP) + PRIMELISTSIZE1 - 1;
101 largest = prime_list[index] + 2;
102 offset -= (index-PRIMELISTSIZE1+1)*STEP + PRIMELISTSIZE1 - 1;
103 }
104
105
106 while(need > obtained && largest > MINPRIME)
107 {
108 isprime = YES;
109 largest -= 2;
110 for(i=0; i<num_prime; i++)
111 if(largest%primes[i] == 0)
112 {
113 isprime = NO;
114 break;
115 }
116
117 if(isprime == YES && offset > 0)
118 offset--;
119 else if(isprime == YES)
120 prime_array[obtained++] = largest;
121 }
122
123 if(need > obtained)
124 fprintf(stderr,"ERROR: Insufficient number of primes: needed %d, obtained %d\n", need, obtained);
125
126 return obtained;
127 }
128
129
130 #if 0
131 main()
132 {
133 int newprimes[1500], np, i;
134
135 np = getprime(2,newprimes,0);
136 np = getprime(2,newprimes+2,9);
137 np = getprime(2,newprimes+4,12);
138
139 for(i=0; i<6; i++)
140 printf("%d. %d \n", i, newprimes[i]);
141
142 /*while(np--)
143 printf("New primes: %d\n", newprimes[np]);
144
145 np = getprime(5,newprimes);
146
147 printf("%d new primes obtained ...\n", np);
148
149 while(np--)
150 printf("New primes: %d\n", newprimes[np]);*/
151 }
152 #endif
153