File Coverage

random_prime.c
Criterion Covered Total %
statement 64 83 77.1
branch 47 58 81.0
condition n/a
subroutine n/a
pod n/a
total 111 141 78.7


line stmt bran cond sub pod time code
1             #include
2             #include
3             #include "csprng.h"
4             #include "primality.h"
5             #include "util.h"
6             #include "prime_nth_count.h"
7             #include "lmo.h"
8             #include "mulmod.h"
9             #include "constants.h"
10             #include "random_prime.h"
11              
12 260           UV random_nbit_prime(void* ctx, UV b)
13             {
14 260           uint32_t start = 0, range;
15             UV n, p;
16 260           switch (b) {
17             case 0:
18 0           case 1: return 0;
19 4 100         case 2: return urandomb(ctx,1) ? 2 : 3;
20 4 100         case 3: return urandomb(ctx,1) ? 5 : 7;
21 4 100         case 4: return urandomb(ctx,1) ? 11 : 13;
22 4           case 5: start = 7; range = 5; break;
23 4           case 6: start = 12; range = 7; break;
24 4           case 7: start = 19; range = 13; break;
25 4           case 8: start = 32; range = 23; break;
26 4           case 9: start = 55; range = 43; break;
27 228           default: break;
28             }
29              
30 248 100         if (start)
31 20           return nth_prime(start + urandomm32(ctx,range));
32              
33 228 50         if (b > BITS_PER_WORD)
34 0           return 0;
35              
36             /* Trivial method */
37 228           p = (UVCONST(1) << (b-1)) + 1;
38             while (1) {
39 2809           n = p + (urandomb(ctx,b-2) << 1);
40 2809 100         if (is_prob_prime(n))
41 228           return n;
42 2581           }
43             }
44              
45 22           UV random_ndigit_prime(void* ctx, UV d)
46             {
47             UV lo, hi;
48 22 50         if ( (d == 0) || (BITS_PER_WORD == 32 && d >= 10) || (BITS_PER_WORD == 64 && d >= 20) ) return 0;
    100          
49 19 100         if (d == 1) return nth_prime(1 + urandomm32(ctx,4));
50 18 100         if (d == 2) return nth_prime(5 + urandomm32(ctx,21));
51 17           lo = powmod(10,d-1,UV_MAX)+1;
52 17           hi = 10*lo-11;
53             while (1) {
54 204           UV n = (lo + urandomm64(ctx,hi-lo+1)) | 1;
55 204 100         if (is_prob_prime(n))
56 17           return n;
57 187           }
58             }
59              
60 22017           UV random_prime(void* ctx, UV lo, UV hi)
61             {
62             UV n, oddrange;
63 22017 100         if (lo > hi) return 0;
64             /* Pull edges in to nearest primes */
65 22015 100         lo = (lo <= 2) ? 2 : next_prime(lo-1);
66 22015 50         hi = (hi >= MPU_MAX_PRIME) ? MPU_MAX_PRIME : prev_prime(hi+1);
67 22015 100         if (lo > hi) return 0;
68             /* There must be at least one prime in the range */
69 22011 100         if (!(lo&1)) lo--; /* treat 2 as 1 */
70 22011           oddrange = ((hi-lo)>>1) + 1; /* look for odds */
71             while (1) {
72 130473           n = lo + 2 * urandomm64(ctx, oddrange);
73 130473 100         if (n == 1 || is_prob_prime(n))
    100          
74 22011 100         return (n == 1) ? 2 : n;
75 108462           }
76             }
77              
78             /* Note that 7 chosen bases or the first 12 prime bases are enough
79             * to guarantee sucess. We could choose to limit to those. */
80 5           int is_mr_random(void* ctx, UV n, UV k) {
81 5 50         if (k >= 3*(n/4))
82 0           return is_prob_prime(n);
83              
84             /* TODO: do 16 at a time */
85 23 100         while (k--) {
86 20           UV base = 2 + urandomm64(ctx, n-2);
87 20 100         if (!miller_rabin(n, &base, 1))
88 20           return 0;
89             }
90 3           return 1;
91             }
92              
93 2           UV random_semiprime(void* ctx, UV b) { /* Even split of bits */
94             static const uint16_t small_semi[] = {35,35,49,65,77,91,143,143,169,299,319,341,377,403};
95             UV min, max, n, L, N;
96              
97 2 50         if (b < 4 || b > BITS_PER_WORD)
    50          
98 0           return 0;
99              
100 2           switch (b) {
101 1           case 4: return 9;
102 0           case 5: return 21;
103 0           case 6: return small_semi[ 0 + urandomm32(ctx,3) ];
104 0           case 7: return small_semi[ 3 + urandomm32(ctx,3) ];
105 0           case 8: return small_semi[ 6 + urandomm32(ctx,3) ];
106 0           case 9: return small_semi[ 9 + urandomm32(ctx,5) ];
107 1           default: break;
108             }
109              
110 1           min = UVCONST(1) << (b-1);
111 1           max = min + (min-1);
112 1           L = b / 2;
113 1           N = b - L;
114              
115             do {
116 4           n = random_nbit_prime(ctx,L) * random_nbit_prime(ctx,N);
117 4 100         } while (n < min || n > max);
    50          
118 1           return n;
119             }
120              
121 1           UV random_unrestricted_semiprime(void* ctx, UV b) { /* generic semiprime */
122             static const unsigned char small_semi[] = {4,6,9,10,14,15,21,22,25,26,33,34,35,38,39,46,49,51,55,57,58,62,65,69,74,77,82,85,86,87,91,93,94,95,106,111,115,118,119,121,122,123};
123             UV min, n;
124              
125 1 50         if (b < 3 || b > BITS_PER_WORD)
    50          
126 0           return 0;
127              
128 1           switch (b) {
129 1           case 3: return small_semi[ 0 + urandomm32(ctx, 2) ];
130 0           case 4: return small_semi[ 2 + urandomm32(ctx, 4) ];
131 0           case 5: return small_semi[ 6 + urandomm32(ctx, 4) ];
132 0           case 6: return small_semi[ 10 + urandomm32(ctx,12) ];
133 0           case 7: return small_semi[ 22 + urandomm32(ctx,20) ];
134 0           default: break;
135             }
136             /* There are faster ways to generate if we could be lax on distribution.
137             * Picking a random prime followed by a second that makes a semiprime in
138             * the range seems obvious and is fast, but the distribution is wrong.
139             * With that method, some semiprimes are much more likely than others. */
140 0           min = UVCONST(1) << (b-1);
141             do {
142 0           n = min + urandomb(ctx,b-1);
143 0 0         } while (!is_semiprime(n));
144 0           return n;
145             }