1
/*
2
*******************************************************************************
3
\file pri.c
4
\brief Prime numbers
5
\project bee2 [cryptographic library]
6
\author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]
7
\created 2012.08.13
8
\version 2016.07.15
9
\license This program is released under the GNU General Public License 
10
version 3. See Copyright Notices in bee2/info.h.
11
*******************************************************************************
12
*/
13

14
#include "bee2/core/mem.h"
15
#include "bee2/core/prng.h"
16
#include "bee2/core/util.h"
17
#include "bee2/core/word.h"
18
#include "bee2/math/pri.h"
19
#include "bee2/math/ww.h"
20
#include "bee2/math/zm.h"
21

22
/*
23
*******************************************************************************
24
Факторная база: первые 1024 нечетных простых
25

26
Построены в Mathematica: 
27
\code
28
	n = 1024;
29
	Table[Prime[i], {i, 2, n + 1}]
30
\endcode
31
*******************************************************************************
32
*/
33

34
static const word _base[] =
35
{
36
	3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 
37
	73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 
38
	157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 
39
	239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 
40
	331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 
41
	421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 
42
	509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 
43
	613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 
44
	709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 
45
	821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 
46
	919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 
47
	1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 
48
	1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 
49
	1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 
50
	1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 
51
	1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 
52
	1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 
53
	1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 
54
	1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 
55
	1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 
56
	1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 
57
	1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 
58
	1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 
59
	2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 
60
	2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 
61
	2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 
62
	2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 
63
	2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 
64
	2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 
65
	2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 
66
	2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 
67
	2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 
68
	2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 
69
	3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 
70
	3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 
71
	3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 
72
	3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 
73
	3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 
74
	3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 
75
	3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 
76
	3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 
77
	3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 
78
	3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 
79
	4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 
80
	4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 
81
	4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 
82
	4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 
83
	4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 
84
	4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 
85
	4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 
86
	4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 
87
	4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 
88
	4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 
89
	5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 
90
	5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 
91
	5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 
92
	5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 
93
	5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 
94
	5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 
95
	5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 
96
	5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 
97
	5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923,	5927, 5939, 5953, 5981, 
98
	5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 
99
	6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 
100
	6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 
101
	6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 
102
	6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 
103
	6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 
104
	6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 
105
	6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 
106
	6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 
107
	6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 
108
	7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 
109
	7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 
110
	7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 
111
	7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 
112
	7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 
113
	7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 
114
	7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 
115
	7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 
116
	7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 
117
	8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 
118
	8147, 8161, 8167, 
119
};
120

121 1
size_t priBaseSize()
122
{
123 1
	return COUNT_OF(_base);
124
}
125

126 1
word priBasePrime(size_t i)
127
{
128 1
	ASSERT(i < priBaseSize());
129 1
	return _base[i];
130
}
131

132
/*
133
*******************************************************************************
134
Произведения последовательных простых из факторной базы: множители 
135
накапливаются до тех пор, пока произведение помещается в машинное слово.
136
Тривиальные произведения из одного множителя не фиксируются.
137

138
Построены в Mathematica: 
139
\code
140
	n = 1024;
141
	w = B_PER_W;
142
	prods = {};
143
	For[i = 2, i <= n + 1, 
144
		For[p = 1; l = 0, 
145
			i + l <= n + 1 && p * Prime[i + ++l] < 2^w, 
146
			p *= Prime[i + l]];
147
		i += l; 
148
		If[l > 1, AppendTo[prods, {p, l}]
149
	]
150
	prods
151
\endcode
152
*******************************************************************************
153
*/
154

155
typedef struct pri_prod_t
156
{
157
	word prod;		/*< произведение простых из факторной базы */
158
	size_t num;		/*< количество множителей */	
159
} pri_prod_t;
160

161
static const pri_prod_t _prods[] =
162
{
163
#if (B_PER_W == 16)
164
	{15015u, 5}, {7429u, 3}, {33263u, 3}, {1763u, 2}, {2491u, 2}, 
165
	{3599u, 2}, {4757u, 2}, {5767u, 2}, {7387u, 2}, {9797u, 2}, {11021u, 2}, 
166
	{12317u, 2}, {16637u, 2}, {19043u, 2}, {22499u, 2}, {25591u, 2}, 
167
	{28891u, 2}, {32399u, 2}, {36863u, 2}, {39203u, 2}, {47053u, 2}, 
168
	{51983u, 2}, {55687u, 2}, {60491u, 2},
169
#elif (B_PER_W == 32)
170
	{3234846615u, 9}, {95041567u, 5}, {907383479u, 5}, {4132280413u, 5},
171
	{121330189u, 4}, {257557397u, 4}, {490995677u, 4}, {842952707u, 4},
172
	{1314423991u, 4}, {2125525169u, 4}, {3073309843u, 4}, {16965341u, 3},
173
	{20193023u, 3}, {23300239u, 3}, {29884301u, 3}, {35360399u, 3},
174
	{42749359u, 3}, {49143869u, 3}, {56466073u, 3}, {65111573u, 3},
175
	{76027969u, 3}, {84208541u, 3}, {94593973u, 3}, {103569859u, 3},
176
	{119319383u, 3}, {133390067u, 3}, {154769821u, 3}, {178433279u, 3},
177
	{193397129u, 3}, {213479407u, 3}, {229580147u, 3}, {250367549u, 3},
178
	{271661713u, 3}, {293158127u, 3}, {319512181u, 3}, {357349471u, 3},
179
	{393806449u, 3}, {422400701u, 3}, {452366557u, 3}, {507436351u, 3},
180
	{547978913u, 3}, {575204137u, 3}, {627947039u, 3}, {666785731u, 3},
181
	{710381447u, 3}, {777767161u, 3}, {834985999u, 3}, {894826021u, 3},
182
	{951747481u, 3}, {1019050649u, 3}, {1072651369u, 3}, {1125878063u, 3},
183
	{1185362993u, 3}, {1267745273u, 3}, {1322520163u, 3}, {1391119619u, 3},
184
	{1498299287u, 3}, {1608372013u, 3}, {1700725291u, 3}, {1805418283u, 3},
185
	{1871456063u, 3}, {2008071007u, 3}, {2115193573u, 3}, {2178429527u, 3},
186
	{2246284699u, 3}, {2385788087u, 3}, {2591986471u, 3}, {2805004793u, 3},
187
	{2922149239u, 3}, {3021320083u, 3}, {3118412617u, 3}, {3265932301u, 3},
188
	{3332392423u, 3}, {3523218343u, 3}, {3711836171u, 3}, {3837879163u, 3},
189
	{3991792529u, 3}, {4139646463u, 3}, {4233155587u, 3}, {2663399u, 2},
190
	{2755591u, 2}, {2782223u, 2}, {2873021u, 2}, {2903591u, 2}, {2965283u, 2},
191
	{3017153u, 2}, {3062491u, 2}, {3125743u, 2}, {3186221u, 2}, {3221989u, 2},
192
	{3301453u, 2}, {3381857u, 2}, {3474487u, 2}, {3504383u, 2}, {3526883u, 2},
193
	{3590989u, 2}, {3648091u, 2}, {3732623u, 2}, {3802499u, 2}, {3904567u, 2},
194
	{3960091u, 2}, {3992003u, 2}, {4028033u, 2}, {4088459u, 2}, {4137131u, 2},
195
	{4235339u, 2}, {4305589u, 2}, {4347221u, 2}, {4384811u, 2}, {4460543u, 2},
196
	{4536899u, 2}, {4575317u, 2}, {4613879u, 2}, {4708819u, 2}, {4862021u, 2},
197
	{4915073u, 2}, {5008643u, 2}, {5048993u, 2}, {5143823u, 2}, {5184713u, 2},
198
	{5244091u, 2}, {5303773u, 2}, {5391563u, 2}, {5475599u, 2}, {5517797u, 2},
199
	{5588447u, 2}, {5659637u, 2}, {5692987u, 2}, {5740807u, 2}, {5827387u, 2},
200
	{5904851u, 2}, {5973127u, 2}, {6066353u, 2}, {6125621u, 2}, {6310063u, 2},
201
	{6426209u, 2}, {6482107u, 2}, {6522907u, 2}, {6682189u, 2}, {6765137u, 2},
202
	{6859157u, 2}, {6969551u, 2}, {7064963u, 2}, {7112873u, 2}, {7182391u, 2},
203
	{7225343u, 2}, {7268407u, 2}, {7338677u, 2}, {7376647u, 2}, {7452899u, 2},
204
	{7535009u, 2}, {7617551u, 2}, {7745053u, 2}, {7806427u, 2}, {7851203u, 2},
205
	{7986227u, 2}, {8065591u, 2}, {8145307u, 2}, {8236819u, 2}, {8363639u, 2},
206
	{8444827u, 2}, {8538059u, 2}, {8678867u, 2}, {8761591u, 2}, {8820899u, 2},
207
	{8999999u, 2}, {9090209u, 2}, {9180851u, 2}, {9272009u, 2}, {9388087u, 2},
208
	{9492557u, 2}, {9603701u, 2}, {9734399u, 2}, {9922331u, 2}, {10036223u, 2}, 
209
	{10137847u, 2}, {10220773u, 2}, {10323353u, 2}, {10400609u, 2}, 
210
	{10575503u, 2}, {10614563u, 2}, {10791029u, 2}, {10916407u, 2}, 
211
	{10995847u, 2}, {11062267u, 2}, {11135533u, 2}, {11242573u, 2}, 
212
	{11329931u, 2}, {11431097u, 2}, {11553137u, 2}, {11716829u, 2}, 
213
	{11923193u, 2}, {11985443u, 2}, {12027023u, 2}, {12215009u, 2},
214
	{12348187u, 2}, {12446783u, 2}, {12503287u, 2}, {12559927u, 2}, 
215
	{12659363u, 2}, {12787751u, 2}, {12873719u, 2}, {13032091u, 2}, 
216
	{13104391u, 2}, {13205947u, 2}, {13329737u, 2}, {13483583u, 2}, 
217
	{13571807u, 2}, {13682597u, 2}, {13793771u, 2}, {13912891u, 2}, 
218
	{14062379u, 2}, {14197823u, 2}, {14333747u, 2}, {14439991u, 2},
219
	{14607683u, 2}, {14745551u, 2}, {14837903u, 2}, {14976851u, 2}, 
220
	{15093209u, 2}, {15280277u, 2}, {15350723u, 2}, {15413467u, 2}, 
221
	{15499933u, 2}, {15657749u, 2}, {15959989u, 2}, {16040021u, 2}, 
222
	{16128247u, 2}, {16192567u, 2}, {16402499u, 2}, {16524161u, 2}, 
223
	{16687189u, 2}, {16777207u, 2}, {16966097u, 2}, {17065157u, 2},
224
	{17189267u, 2}, {17288963u, 2}, {17547577u, 2}, {17757787u, 2}, 
225
	{17842151u, 2}, {17943671u, 2}, {18045479u, 2}, {18147599u, 2}, 
226
	{18249983u, 2}, {18369787u, 2}, {18593119u, 2}, {18818243u, 2}, 
227
	{18948593u, 2}, {19079399u, 2}, {19307227u, 2}, {19492189u, 2}, 
228
	{19642543u, 2}, {19793597u, 2}, {19891591u, 2}, {20088323u, 2},
229
	{20249951u, 2}, {20385221u, 2}, {20439437u, 2}, {20684303u, 2}, 
230
	{20830087u, 2}, {21040553u, 2}, {21159991u, 2}, {21427577u, 2}, 
231
	{21538877u, 2}, {21622499u, 2}, {21715591u, 2}, {21864967u, 2}, 
232
	{22061773u, 2}, {22297283u, 2}, {22382357u, 2}, {22610009u, 2}, 
233
	{22896221u, 2}, {22953677u, 2}, {23039999u, 2}, {23184221u, 2},
234
	{23483491u, 2}, {23755867u, 2}, {23970767u, 2}, {24147371u, 2}, 
235
	{24324623u, 2}, {24403591u, 2}, {24542107u, 2}, {24681023u, 2}, 
236
	{24800351u, 2}, {24960007u, 2}, {25060027u, 2}, {25160231u, 2}, 
237
	{25310897u, 2}, {25553009u, 2}, {25796237u, 2}, {25938613u, 2}, 
238
	{26050807u, 2}, {26173447u, 2}, {26522491u, 2}, {26718557u, 2},
239
	{26873831u, 2}, {27071173u, 2}, {27342437u, 2}, {27405221u, 2}, 
240
	{27741253u, 2}, {27878399u, 2}, {28089991u, 2}, {28259807u, 2}, 
241
	{28515551u, 2}, {28793731u, 2}, {29052091u, 2}, {29192393u, 2}, 
242
	{29322221u, 2}, {29430589u, 2}, {29582717u, 2}, {29658907u, 2}, 
243
	{29964667u, 2}, {30041357u, 2}, {30272003u, 2}, {30393133u, 2},
244
	{30514567u, 2}, {30735767u, 2}, {30980347u, 2}, {31102913u, 2}, 
245
	{31438193u, 2}, {31809599u, 2}, {31911197u, 2}, {31979021u, 2}, 
246
	{32080871u, 2}, {32330587u, 2}, {32455793u, 2}, {32649787u, 2}, 
247
	{32936117u, 2}, {33016507u, 2}, {33419957u, 2}, {33593591u, 2}, 
248
	{33756091u, 2}, {33918967u, 2}, {34117277u, 2}, {34222499u, 2},
249
	{34327877u, 2}, {34433423u, 2}, {34574399u, 2}, {34809991u, 2}, 
250
	{35105621u, 2}, {35354867u, 2}, {35808247u, 2}, {36108077u, 2}, 
251
	{36397073u, 2}, {36542021u, 2}, {36723551u, 2}, {36917767u, 2}, 
252
	{37088099u, 2}, {37295413u, 2}, {37527851u, 2}, {37675019u, 2}, 
253
	{37908613u, 2}, {38254081u, 2}, {38452397u, 2}, {38613787u, 2},
254
	{38750609u, 2}, {39087479u, 2}, {39262747u, 2}, {39363067u, 2}, 
255
	{39601813u, 2}, {39765611u, 2}, {39942391u, 2}, {40106873u, 2}, 
256
	{40297079u, 2}, {40449599u, 2}, {40576891u, 2}, {40755431u, 2}, 
257
	{41075137u, 2}, {41447723u, 2}, {41731519u, 2}, {41951513u, 2}, 
258
	{42327811u, 2}, {42745363u, 2}, {42928703u, 2}, {43112347u, 2},
259
	{43217467u, 2}, {43428019u, 2}, {43731733u, 2}, {44155961u, 2}, 
260
	{44355599u, 2}, {44568967u, 2}, {44756099u, 2}, {44916803u, 2}, 
261
	{45077771u, 2}, {45360221u, 2}, {45724643u, 2}, {45968399u, 2}, 
262
	{46131263u, 2}, {46416869u, 2}, {46621583u, 2}, {46744553u, 2}, 
263
	{47059591u, 2}, {47196899u, 2}, {47485817u, 2}, {47734277u, 2},
264
	{48052399u, 2}, {48358091u, 2}, {48497287u, 2}, {48636667u, 2}, 
265
	{48818153u, 2}, {48985997u, 2}, {49224247u, 2}, {49463053u, 2}, 
266
	{49702451u, 2}, {50041451u, 2}, {50495227u, 2}, {50751367u, 2}, 
267
	{50979479u, 2}, {51380143u, 2}, {51696091u, 2}, {51969677u, 2}, 
268
	{52070647u, 2}, {52316273u, 2}, {52490021u, 2}, {52823599u, 2},
269
	{53319179u, 2}, {53509189u, 2}, {53758223u, 2}, {54022499u, 2}, 
270
	{54479017u, 2}, {54967387u, 2}, {55383283u, 2}, {55621763u, 2}, 
271
	{55935437u, 2}, {56070143u, 2}, {56294993u, 2}, {56550391u, 2}, 
272
	{56746073u, 2}, {56911927u, 2}, {57062891u, 2}, {57259453u, 2}, 
273
	{57456391u, 2}, {57608099u, 2}, {57836021u, 2}, {58216819u, 2},
274
	{58461307u, 2}, {58844237u, 2}, {59043847u, 2}, {59213009u, 2}, 
275
	{59444051u, 2}, {59675621u, 2}, {60015973u, 2}, {60186563u, 2}, 
276
	{60699677u, 2}, {61152391u, 2}, {61387189u, 2}, {61779551u, 2}, 
277
	{62015621u, 2}, {62110157u, 2}, {62473207u, 2}, {62773913u, 2}, 
278
	{62964221u, 2}, {63202499u, 2}, {63648259u, 2}, {64160099u, 2},
279
	{64448663u, 2}, {64899127u, 2}, {65205589u, 2}, {65415743u, 2}, 
280
	{65561393u, 2}, {65836987u, 2}, {66178081u, 2}, {66650887u, 2},
281
#else
282
	{16294579238595022365u, 15}, {7145393598349078859u, 10},
283
	{6408001374760705163u, 9}, {690862709424854779u, 8},
284
	{4312024209383942993u, 8}, {71235931512604841u, 7},
285
	{192878245514479103u, 7}, {542676746453092519u, 7},
286
	{1230544604996048471u, 7}, {2618501576975440661u, 7},
287
	{4771180125133726009u, 7}, {9247077179230889629u, 7},
288
	{32156968791364271u, 6}, {46627620659631719u, 6}, {64265583549260393u, 6}, 
289
	{88516552714582021u, 6}, {131585967012906751u, 6}, 
290
	{182675399263485151u, 6}, {261171077386532413u, 6}, 
291
	{346060227726080771u, 6}, {448604664249794309u, 6},
292
	{621993868801161359u, 6}, {813835565706097817u, 6}, 
293
	{1050677302683430441u, 6}, {1294398862104002783u, 6}, 
294
	{1615816556891330179u, 6}, {1993926996710486603u, 6}, 
295
	{2626074105497143999u, 6}, {3280430033433832817u, 6}, 
296
	{4076110663011485663u, 6}, {4782075577404875363u, 6}, 
297
	{5906302864496324923u, 6}, {7899206880638488339u, 6}, 
298
	{9178333502078117453u, 6}, {10680076322389870367u, 6}, 
299
	{12622882367374918799u, 6}, {14897925470078818423u, 6}, 
300
	{17264316336968551717u, 6}, {11896905306684389u, 5}, 
301
	{13580761294555417u, 5}, {15289931661301991u, 5}, 
302
	{17067874133764579u, 5}, {19008757261780379u, 5}, 
303
	{21984658219193689u, 5}, {23721541361298551u, 5}, 
304
	{26539432378378657u, 5}, {30167221680049747u, 5}, 
305
	{32433198277139683u, 5}, {35517402656173043u, 5}, 
306
	{39100537712055041u, 5}, {42477532426853543u, 5}, 
307
	{45618621452253523u, 5}, {52071972962579407u, 5}, 
308
	{57329264013213233u, 5}, {61692083285823527u, 5}, 
309
	{66885169838978461u, 5}, {72186879569637319u, 5}, 
310
	{77103033998665567u, 5}, {82549234838454463u, 5},
311
	{89609394623390063u, 5}, {100441814079170659u, 5}, 
312
	{109045745121501371u, 5}, {120230527473437819u, 5}, 
313
	{131125107904515419u, 5}, {138612182127286823u, 5}, 
314
	{144712752835963307u, 5}, {152692680370726429u, 5}, 
315
	{164664356404541573u, 5}, {175376065798883557u, 5}, 
316
	{187958301132741257u, 5}, {203342285718459187u, 5}, 
317
	{219115706321995421u, 5}, {235226887496676263u, 5}, 
318
	{253789253193479219u, 5}, {271717583502831491u, 5}, 
319
	{293266389497362763u, 5}, {321821627692439603u, 5}, 
320
	{339856237957830049u, 5}, {362469273063260281u, 5}, 
321
	{390268963330916339u, 5}, {408848490015359209u, 5}, 
322
	{429644565036857699u, 5}, {458755816747679897u, 5}, 
323
	{495450768525623033u, 5}, {523240424009891327u, 5}, 
324
	{551070603968128061u, 5}, {574205321266688311u, 5}, 
325
	{606829434176923693u, 5}, {637763212653336997u, 5}, 
326
	{676538378976146257u, 5}, {710263471119657661u, 5}, 
327
	{754496879875465343u, 5}, {800075738315885429u, 5}, 
328
	{845197573085733239u, 5}, {894146362391888161u, 5}, 
329
	{930105507041885771u, 5}, {985345849616172623u, 5}, 
330
	{1040222328124784927u, 5}, {1091468150538871153u, 5}, 
331
	{1150933747479716653u, 5}, {1210604027868555713u, 5}, 
332
	{1277530693373553361u, 5}, {1350088100087645657u, 5}, 
333
	{1398676120233167591u, 5}, {1459450139327525269u, 5}, 
334
	{1555755169940697937u, 5}, {1645735334920325819u, 5}, 
335
	{1732866357938791147u, 5}, {1815492864312158099u, 5}, 
336
	{1894564116319543619u, 5}, {1993720023757886939u, 5}, 
337
	{2103356633892712673u, 5}, {2180099035358103487u, 5}, 
338
	{2277315690161244011u, 5}, {2390146379836558999u, 5}, 
339
	{2522126262040806983u, 5}, {2613887383383648311u, 5}, 
340
	{2795412145600606001u, 5}, {2919958494381348367u, 5}, 
341
	{3012266980379247553u, 5}, {3119360766859522543u, 5}, 
342
	{3216618305468232557u, 5}, {3385071039891962579u, 5}, 
343
	{3509427475807939163u, 5}, {3700008514672760651u, 5}, 
344
	{3873423910591589033u, 5}, {4050200067084600439u, 5}, 
345
	{4233429923647833421u, 5}, {4472862244562412787u, 5}, 
346
	{4638587045132438407u, 5}, {4765110805097342489u, 5}, 
347
	{4951886102290887619u, 5}, {5103665412856065733u, 5}, 
348
	{5306636943410213377u, 5}, {5581202660702121667u, 5}, 
349
	{5774946339890457283u, 5}, {5948565823654343479u, 5}, 
350
	{6175776426345604697u, 5}, {6454381412132929663u, 5}, 
351
	{6685489970462824483u, 5}, {6864273057227912189u, 5}, 
352
	{7020466399135670969u, 5}, {7326535375987923521u, 5}, 
353
	{7795297680723533551u, 5}, {8101368883577379127u, 5}, 
354
	{8353548973662446233u, 5}, {8642941459163335097u, 5}, 
355
	{8989536585105548947u, 5}, {9281597047973093449u, 5}, 
356
	{9623989560555822323u, 5}, {9884958654427267811u, 5}, 
357
	{10161260209201654649u, 5}, {10427299413974952277u, 5}, 
358
	{10759022378261015069u, 5}, {11290266633494854903u, 5}, 
359
	{11852839900193264543u, 5}, {12209591760462366097u, 5}, 
360
	{12604874462407914499u, 5}, {13152204116524238149u, 5}, 
361
	{13487108616207513373u, 5}, {13935742926215786057u, 5}, 
362
	{14426307807825845411u, 5}, {14869398407236523377u, 5}, 
363
	{15287602532539390847u, 5}, {15824557271701228177u, 5}, 
364
	{16348641675671844607u, 5}, {16684838939207505557u, 5}, 
365
	{17148166373867218913u, 5}, {17832011695956938489u, 5}, 
366
	{2587278248197793u, 4}, {2656152548121013u, 4}, {2706094705771019u, 4}, 
367
	{2746082268411733u, 4}, {2816510930505221u, 4}, {2876558914811147u, 4}, 
368
	{2943092641403483u, 4}, {3044274349991521u, 4}, {3111227620115431u, 4}, 
369
	{3156468307693999u, 4}, {3209012615864543u, 4}, {3247559087000957u, 4}, 
370
	{3289921520014123u, 4}, {3331823223534079u, 4}, {3403431328122433u, 4}, 
371
	{3474390126259739u, 4}, {3519861126859459u, 4}, {3581490458694233u, 4}, 
372
	{3653304933840151u, 4}, {3753973384118899u, 4}, {3831297220366171u, 4}, 
373
	{3880220695063499u, 4}, {3952510531166773u, 4}, {4022729025799241u, 4}, 
374
	{4135032598497637u, 4}, {4231785801620803u, 4}, {4288747235209999u, 4}, 
375
	{4356965458481947u, 4}, {66650887u, 2},
376
#endif
377
};
378

379 1
void priBaseMod(word mods[], const word a[], size_t n, size_t count)
380
{
381
	size_t i, j;
382
	// pre
383 1
	ASSERT(wwIsValid(a, n));
384 1
	ASSERT(count <= priBaseSize());
385 1
	ASSERT(wwIsValid(mods, count));
386
	// пробегаем произведения простых из факторной базы
387 1
	for (i = j = 0; i < count && j < COUNT_OF(_prods); ++j)
388
	{
389 1
		size_t num = _prods[j].num;
390
		// t <- a mod _base[i] * ... * _base[i + num - 1]
391 1
		word t = zzModW(a, n, _prods[j].prod);
392
		// mods[i] <- t mod _base[i]
393 1
		while (num-- && i < count)
394 1
			mods[i] = t % _base[i], ++i;
395
	}
396
	// пробегаем оставшиеся простые из факторной базы
397 1
	for (; i < count; ++i)
398 0
		mods[i] = zzModW(a, n, _base[i]);
399
}
400

401
/*
402
*******************************************************************************
403
Использование факторной базы
404

405
\todo Алгоритм Бернштейна выделения гладкой части 
406
[http://cr.yp.to/factorization/smoothparts-20040510.pdf]:
407
	z <- p1 p2 .... pm \mod a
408
	y <- z^{2^e} \mod a (e --- минимальное натуральное т.ч. 3^{2^e} > a)
409
	return \gcd(x, y)
410
*******************************************************************************
411
*/
412

413 1
bool_t priIsSieved(const word a[], size_t n, size_t base_count, void* stack)
414
{
415
	// переменные в stack
416
	word* mods;
417
	// pre
418 1
	ASSERT(base_count <= priBaseSize());
419
	// четное?
420 1
	n = wwWordSize(a, n);
421 1
	if (zzIsEven(a, n))
422 0
		return FALSE;
423
	// малое a?
424 1
	if (n == 1)
425
		// при необходимости скоррректировать факторную базу
426 0
		while (base_count > 0 && priBasePrime(base_count - 1) > a[0])
427 0
			--base_count;
428
	// раскладка stack
429 1
	mods = (word*)stack;
430 1
	stack = mods + base_count;
431
	// найти остатки
432 1
	priBaseMod(mods, a, n, base_count);
433
	// есть нулевые остатки?
434 1
	while (base_count--)
435 1
		if (mods[base_count] == 0)
436 1
			return FALSE;
437
	// нет
438 1
	return TRUE;
439
}
440

441 1
size_t priIsSieved_deep(size_t base_count)
442
{
443 1
	return O_OF_W(base_count);
444
}
445

446 1
bool_t priIsSmooth(const word a[], size_t n, size_t base_count, void* stack)
447
{
448
	register size_t i;
449
	register word mod;
450
	// переменные в stack
451 1
	word* t = (word*)stack;
452 1
	stack = t + n;
453
	// pre
454 1
	ASSERT(base_count <= priBaseSize());
455
	// t <- a 
456 1
	wwCopy(t, a, n);
457
	// разделить t на степень 2
458 1
	i = wwLoZeroBits(t, n);
459 1
	wwShLo(t, n, i);
460 1
	n = wwWordSize(t, n);
461 1
	if (wwIsZero(t, n))
462
	{
463 0
		i = 0;
464 0
		return TRUE;
465
	}
466
	// цикл по простым из факторной базы
467 1
	for (i = 0; i < base_count;)
468
	{
469 1
		mod = _base[i] < WORD_BIT_HALF ? 
470 1
			zzModW2(t, n, _base[i]) : zzModW(t, n, _base[i]);
471
		// делится на простое?
472 1
		if (mod == 0)
473
		{
474 1
			zzDivW(t, t, n, _base[i]);
475 1
			n = wwWordSize(t, n);
476 1
			if (wwIsZero(t, n))
477
			{
478 0
				i = 0, mod = 0;
479 0
				return TRUE;
480
			}
481
		}
482
		// .. не делится
483
		else 
484 1
			++i;
485
	}
486 1
	i = 0, mod = 0;
487 1
	return FALSE;
488
}
489

490 1
size_t priIsSmooth_deep(size_t n)
491
{
492 1
	return O_OF_W(n);
493
}
494

495
/*
496
*******************************************************************************
497
Проверка простоты малых чисел
498

499
Применяется тест Миллера --- Рабина со специально подобранными основаниями.
500
Успешное завершение теста на всех основаниях из списка _base16 гарантирует
501
простоту чисел вплоть до 1373653, из списка _base32 --- вплоть
502
до 4759123141, из списка _base64 --- для всех 64-разрядных чисел.
503

504
Список оснований: http://miller-rabin.appspot.com (список обновляется).
505

506
\todo Проверка на малые делители позволит уменьшить число оснований / 
507
сократить время проверки?
508
*******************************************************************************
509
*/
510

511
const word _bases16[] = {2, 3};
512
const word _bases32[] = {2, 7, 61};
513
#if (B_PER_W == 64)
514
const word _bases64[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
515
#endif
516

517 1
bool_t priIsPrimeW(register word a, void* stack)
518
{
519
	const word* bases;
520
	register word r;
521
	register size_t s;
522
	register size_t iter;
523
	register size_t i;
524
	register word base;
525
	register dword prod;
526
	// маленькое или четное?
527 1
	if (a <= 3 || a % 2 == 0)
528 1
		return a == 2 || a == 3;
529
	// a - 1 = 2^s r (r -- нечетное)
530 1
	for (r = a - 1, s = 0; r % 2 == 0; r >>= 1, ++s);
531 1
	ASSERT(s > 0 && WORD_BIT_POS(s) * r + 1 == a);
532
	// выбираем базис
533
#if (B_PER_W == 16)
534
	bases = _bases16, iter = COUNT_OF(_bases16);
535
#elif (B_PER_W == 32)
536
	if (a < 1373653)
537
		bases = _bases16, iter = COUNT_OF(_bases16);
538
	else
539
		bases = _bases32, iter = COUNT_OF(_bases32);
540
#elif (B_PER_W == 64)
541 1
	if (a < 1373653)
542 1
		bases = _bases16, iter = COUNT_OF(_bases16);
543 0
	else if (a < 4759123141)
544 0
		bases = _bases32, iter = COUNT_OF(_bases32);
545
	else
546 0
		bases = _bases64, iter = COUNT_OF(_bases64);
547
#endif
548
	// итерации
549 1
	while (iter--)
550
	{
551
		// _bases[iter]^r \equiv \pm 1 \mod a?
552 1
		base = zzPowerModW(bases[iter], r, a, stack);
553 1
		if (base == 1 || base == a - 1)
554 1
			continue;
555
		// base^{2^i} \equiv - 1\mod a?
556 1
		for (i = s - 1; i--;)
557
		{
558 1
			prod = base;
559 1
			prod *= base, base = prod % a;
560 1
			if (base == a - 1)
561 1
				break;
562 1
			if (base == 1)
563
			{
564 0
				r = base = 0, s = iter = i = 0, prod = 0;
565 0
				return FALSE;
566
			}
567
		}
568 1
		if (i == SIZE_MAX)
569
		{
570 1
			r = base = 0, s = iter = i = 0, prod = 0;
571 1
			return FALSE;
572
		}
573
	}
574 1
	r = base = 0, s = iter = i = 0, prod = 0;
575 1
	return TRUE;
576
}
577

578 1
size_t priIsPrimeW_deep()
579
{
580 1
	return zzPowerModW_deep();
581
}
582

583
/*
584
*******************************************************************************
585
Тест Рабина -- Миллера
586

587
При операциях \mod a может использоваться умножение Монтгомери. При этом
588
случайное основание выбирается отличным от R \mod a, -R \mod a (R = B^n),
589
а получаемые степени сравниваются не с \pm 1 \mod a, а с \pm R \mod a.
590

591
В начале функции проверяется, что a >= 49. Поэтому основание base будет 
592
совпадать с \pm R \mod a с вероятностью p < 2/(49 - 1) = 1/24 и число попыток 
593
генерации можно ограничить величиной 
594
B_PER_IMPOSSLIBLE / log_2(24) < B_PER_IMPOSSLIBLE / 4.5.
595
*******************************************************************************
596
*/
597

598 1
bool_t priRMTest(const word a[], size_t n, size_t iter, void* stack)
599
{
600
	register size_t s;
601
	register size_t m;
602
	register size_t i;
603
	// переменные в stack
604 1
	word* r = (word*)stack;
605 1
	word* base = r + n;
606 1
	qr_o* qr = (qr_o*)(base + n);
607 1
	octet* combo_state = (octet*)qr + zmCreate_keep(O_OF_W(n));
608 1
	stack = combo_state + prngCOMBO_keep();
609
	// pre
610 1
	ASSERT(wwIsValid(a, n));
611
	// нормализация
612 1
	n = wwWordSize(a, n);
613
	// четное?
614 1
	if (zzIsEven(a, n))
615 0
		return wwCmpW(a, n, 2) == 0;
616
	// маленькое?
617 1
	if (n == 1 && a[0] < 49)
618 1
		return a[0] != 1 && (a[0] == 3 || a[0] % 3) && (a[0] == 5 || a[0] % 5);
619
	// подготовить генератор
620 1
	prngCOMBOStart(combo_state, utilNonce32());
621
	// создать кольцо
622 1
	wwTo(base, O_OF_W(n), a);
623 1
	zmCreate(qr, (octet*)base, memNonZeroSize(base, O_OF_W(n)), stack);
624
	// a - 1 = r 2^s (r -- нечетное)
625 1
	wwCopy(r, a, n);
626 1
	zzSubW2(r, n, 1);
627 1
	s = wwLoZeroBits(r, n);
628 1
	wwShLo(r, n, s);
629 1
	m = wwWordSize(r, n);
630
	// итерации
631 1
	while (iter--)
632
	{
633
		// base <-R {1, \ldots, a - 1} \ {\pm one}
634 1
		i = 0;
635
		do
636 1
			if (i++ * 45 > B_PER_IMPOSSIBLE * 10 || 
637 1
				!zzRandNZMod(base, a, n, prngCOMBOStepR, combo_state))
638
			{
639 0
				s = m = 0;
640 0
				return FALSE;
641
			}
642 1
		while (wwEq(base, qr->unity, n) || zzIsSumEq(a, base, qr->unity, n));
643
		// base <- base^r \mod a
644 1
		qrPower(base, base, r, m, qr, stack);
645
		// base == \pm one => тест пройден
646 1
		if (wwEq(base, qr->unity, n) ||
647 1
			zzIsSumEq(a, base, qr->unity, n))
648 1
			continue;
649
		// base^{2^i} \equiv -1 \mod a?
650 1
		for (i = s; i--;)
651
		{
652 1
			qrSqr(base, base, qr, stack);
653 1
			if (wwEq(base, qr->unity, n))
654
			{
655 0
				s = m = i = 0;
656 0
				return FALSE;
657
			}
658 1
			if (zzIsSumEq(a, base, qr->unity, n))
659 1
				break;
660
		}
661 1
		if (i == SIZE_MAX)
662
		{
663 1
			s = m = i = 0;
664 1
			return FALSE;
665
		}
666
	}
667 1
	s = m = i = 0;
668
	// простое
669 1
	return TRUE;
670
}
671

672 1
size_t priRMTest_deep(size_t n)
673
{
674 1
	size_t qr_deep = zmCreate_deep(O_OF_W(n));
675 1
	return O_OF_W(2 * n) + zmCreate_keep(O_OF_W(n)) + prngCOMBO_keep() +
676 1
		utilMax(2,
677
			qr_deep,
678
			qrPower_deep(n, n, qr_deep));
679
}
680

681 1
bool_t priIsPrime(const word a[], size_t n, void* stack)
682
{
683 1
	return priRMTest(a, n, (B_PER_IMPOSSIBLE + 1) / 2, stack);
684
}
685

686 1
size_t priIsPrime_deep(size_t n)
687
{
688 1
	return priRMTest_deep(n);
689
}
690

691
/*
692
*******************************************************************************
693
Простое Софи Жермен
694

695
q -- нечетное простое => p = 2q + 1 -- простое <=> [теорема Демитко]
696
1) 2^2 \not\equiv 1 \mod p;
697
2) 2^2q \equiv 1 \mod p.
698
*******************************************************************************
699
*/
700

701 1
bool_t priIsSGPrime(const word q[], size_t n, void* stack)
702
{
703 1
	size_t no = O_OF_W(n + 1);
704
	// переменные в stack
705
	word* p;
706
	qr_o* qr;
707
	// pre
708 1
	ASSERT(zzIsOdd(q, n) && wwCmpW(q, n, 1) > 0);
709
	// раскладка стек
710 1
	p = (word*)stack;
711 1
	qr = (qr_o*)(p + n + 1);
712 1
	stack = (octet*)qr + zmCreate_keep(no);
713
	// p <- 2q + 1
714 1
	wwCopy(p, q, n);
715 1
	p[n] = 0;
716 1
	wwShHi(p, n + 1, 1);
717 1
	++p[0];
718
	// создать кольцо \mod p
719 1
	no = wwOctetSize(p, n + 1);
720 1
	wwTo(p, no, p);
721 1
	zmCreate(qr, (octet*)p, no, stack);
722
	// p <- 4^q (в кольце qr)
723 1
	qrAdd(p, qr->unity, qr->unity, qr);
724 1
	qrAdd(p, p, p, qr);
725 1
	qrPower(p, p, q, n, qr, stack);
726
	// 4^q == 1 (в кольце qr)?
727 1
	return qrCmp(p, qr->unity, qr) == 0;
728
}
729

730 1
size_t priIsSGPrime_deep(size_t n)
731
{
732 1
	const size_t no = O_OF_W(n + 1);
733 1
	const size_t qr_deep = zmCreate_deep(no);
734 1
	return no + zmCreate_keep(no) +
735 1
		utilMax(2,
736
			qr_deep,
737
			qrPower_deep(n + 1, n, qr_deep));
738
}
739

740
/*
741
*******************************************************************************
742
Следующее простое
743
*******************************************************************************
744
*/
745

746 1
bool_t priNextPrimeW(word p[1], register word a, void* stack)
747
{
748
	register size_t l;
749
	// p <- a, l <- битовая длина a
750 1
	p[0] = a, l = wwBitSize(p, 1);
751
	// 0-битовых и 1-битовых простых не существует
752 1
	if (l <= 1)
753 0
		return FALSE;
754
	// сделать p нечетным
755 1
	p[0] |= 1;
756
	// поиск
757 1
	while (!priIsPrimeW(p[0], stack))
758
	{
759 1
		p[0] += 2;
760 1
		if (wwBitSize(p, 1) != l)
761
		{
762 0
			l = 0;
763 0
			return FALSE;
764
		}
765
	}
766 1
	l = 0;
767 1
	return TRUE;
768
}
769

770 1
size_t priNextPrimeW_deep()
771
{
772 1
	return priIsPrimeW_deep();
773
}
774

775 1
bool_t priNextPrime(word p[], const word a[], size_t n, size_t trials,
776
	size_t base_count, size_t iter, void* stack)
777
{
778
	size_t l;
779
	size_t i;
780
	bool_t base_success;
781
	// переменные в stack
782
	word* mods;
783
	// pre
784 1
	ASSERT(wwIsSameOrDisjoint(a, p, n));
785 1
	ASSERT(base_count <= priBaseSize());
786
	// раскладка stack
787 1
	mods = (word*)stack;
788 1
	stack = mods + base_count;
789
	// l <- битовая длина a
790 1
	l = wwBitSize(a, n);
791
	// 0-битовых и 1-битовых простых не существует
792 1
	if (l <= 1)
793 0
		return FALSE;
794
	// p <- минимальное нечетное >= a
795 1
	wwCopy(p, a, n);
796 1
	p[0] |= 1;
797
	// малое p?
798 1
	if (n == 1)
799
		// при необходимости скоррректировать факторную базу
800 1
		while (base_count > 0 && priBasePrime(base_count - 1) >= p[0])
801 0
			--base_count;
802
	// рассчитать остатки от деления на малые простые
803 1
	priBaseMod(mods, p, n, base_count);
804 1
	for (i = 0, base_success = TRUE; i < base_count; ++i)
805 1
		if (mods[i] == 0)
806
		{
807 1
			base_success = FALSE;
808 1
			break;
809
		}
810
	// попытки
811 1
	while (trials == SIZE_MAX || trials--)
812
	{
813
		// проверка простоты
814 1
		if (base_success && priRMTest(p, n, iter, stack))
815 1
			return TRUE;
816
		// к следующему кандидату
817 1
		if (zzAddW2(p, n, 2) || wwBitSize(p, n) > l)
818 1
			return FALSE;
819 1
		for (i = 0, base_success = TRUE; i < base_count; ++i)
820
		{
821 1
			if (mods[i] < _base[i] - 2)
822 1
				mods[i] += 2;
823 1
			else if (mods[i] == _base[i] - 1)
824 1
				mods[i] = 1;
825
			else
826 1
				mods[i] = 0, base_success = FALSE;
827
		}
828
	}
829 0
	return FALSE;
830
}
831

832 1
size_t priNextPrime_deep(size_t n, size_t base_count)
833
{
834 1
	return base_count * O_PER_W + priRMTest_deep(n);
835
}
836

837
/*
838
*******************************************************************************
839
Расширение простого
840

841
Теорема Демитко. Если q -- нечетное простое, p = 2qr + 1, где 2r < 4q + 1,
842
и выполнены условия:
843
1) 2^{2qr} \equiv 1 \mod p;
844
2) 2^{2r} \not\equiv 1 \mod p,
845
то p -- простое.
846

847
\remark Если l <= 2 * lq, где lq = bitlen(q), то условие 2r < 4q + 1 будет
848
выполнено:
849
2r = (p - 1) / q < 2^l / (2^{lq - 1}} = 2^{l - lq + 1} <= 2^{lq + 2} < 4q.
850

851
Построение p:
852
1) t <-R {2^{l - 2} + 1,..., 2^{l - 1} - 1};
853
2) t <- t + 2^{l - 1};
854
3) r <- ceil(t / q);
855
4) p <- 2qr + 1;
856
5) если bitlen(p) != l, то вернуться к шагу 1.
857

858
\remark Если t укладывается в m слов, q -- в n слов, то r на шаге 3)
859
укладывается в m - n + 1 слов. Действительно, максимальное r получается
860
при t = B^m - 1, q = B^{n - 1} и равняется B^{m - n + 1} - 1.
861
*******************************************************************************
862
*/
863

864 1
bool_t priExtendPrime(word p[], size_t l, const word q[], size_t n,
865
	size_t trials, size_t base_count, gen_i rng, void* rng_state, void* stack)
866
{
867 1
	const size_t m = W_OF_B(l);
868 1
	const size_t mo = O_OF_B(l);
869
	size_t i;
870
	// переменные в stack
871
	word* r;
872
	word* t;
873
	word* four;
874
	word* mods;
875
	word* mods1;
876
	qr_o* qr;
877
	// pre
878 1
	ASSERT(wwIsDisjoint2(q, n, p, m));
879 1
	ASSERT(zzIsOdd(q, n) && wwCmpW(q, n, 3) >= 0);
880 1
	ASSERT(wwBitSize(q, n) + 1 <= l && l <= 2 * wwBitSize(q, n));
881 1
	ASSERT(base_count <= priBaseSize());
882 1
	ASSERT(rng != 0);
883
	// подкорректировать n
884 1
	n = wwWordSize(q, n);
885
	// раскладка stack
886 1
	r = (word*)stack;
887 1
	t = r + m - n + 1;
888 1
	four = t + m + 1;
889 1
	mods = four + m;
890 1
	mods1 = mods + base_count;
891 1
	qr = (qr_o*)(mods1 + base_count);
892 1
	stack = (octet*)qr + zmCreate_keep(mo);
893
	// малое p?
894 1
	if (l < B_PER_W)
895
		// при необходимости уменьшить факторную базу
896 1
		while (base_count > 0 && 
897 1
			priBasePrime(base_count - 1) > WORD_BIT_POS(l - 1))
898 0
			--base_count;
899
	// попытки
900 1
	while (trials == SIZE_MAX || trials--)
901
	{
902
		// t <-R [2^{l - 2}, 2^{l - 1})
903 1
		rng(t, mo, rng_state);
904 1
		wwFrom(t, t, mo);
905 1
		wwTrimHi(t, m, l - 2);
906 1
		wwSetBit(t, l - 2, 1);
907
		// r <- t \div q
908 1
		zzDiv(r, t, t, m, q, n, stack);
909 1
		if (!wwIsZero(t, m))
910 1
			VERIFY(zzAddW2(r, m - n + 1, 1) == 0);
911
		// t <- q * r
912 1
		zzMul(t, q, n, r, m - n + 1, stack);
913 1
		if (wwBitSize(t, m) > l - 1)
914 1
			continue;
915
		// p <- 2 * t + 1
916 1
		wwCopy(p, t, m);
917 1
		wwShHi(p, m, 1);
918 1
		++p[0];
919 1
		ASSERT(wwBitSize(p, m) == l);
920
		// рассчитать вычеты p, 2q по малым модулям
921 1
		priBaseMod(mods, p, m, base_count);
922 1
		priBaseMod(mods1, q, n, base_count);
923 1
		for (i = 0; i < base_count; ++i)
924 1
			if ((mods1[i] += mods1[i]) >= _base[i])
925 1
				mods1[i] -= _base[i];
926
		// проверка простоты
927
		while (1)
928
		{
929
			// p делится на малые простые?
930 1
			for (i = 0; i < base_count; ++i)
931 1
				if (mods[i] == 0)
932 1
					break;
933
			// не делится: тест Демитко
934 1
			if (i == base_count)
935
			{
936
				// создать кольцо вычетов \mod p
937 1
				wwTo(t, mo, p);
938 1
				zmCreate(qr, (octet*)t, mo, stack);
939
				// four <- 4 [в кольце qr]
940 1
				qrAdd(four, qr->unity, qr->unity, qr);
941 1
				qrAdd(four, four, four, qr);
942
				// 4^r \mod p != 1?
943 1
				qrPower(t, four, r, m - n + 1, qr, stack);
944 1
				if (qrCmp(t, qr->unity, qr) != 0)
945
				{
946
					// (4^r)^q \mod p == 1?
947 1
					qrPower(t, t, q, n, qr, stack);
948 1
					if (qrCmp(t, qr->unity, qr) == 0)
949 1
						return TRUE;
950
				}
951
			}
952
			// p <- p + 2q, переполнение?
953 1
			if (zzAddW2(p + n, m - n, zzAdd2(p, q, n)) ||
954 1
				zzAddW2(p + n, m - n, zzAdd2(p, q, n)) ||
955 1
				wwBitSize(p, m) > l)
956
				break;
957
			// r <- r + 1, t <- t + q, пересчитать mods
958 1
			zzAddW2(r, m - n + 1, 1);
959 1
			zzAddW2(t + n, m - n, zzAdd2(t, q, n));
960 1
			for (i = 0; i < base_count; ++i)
961 1
				if ((mods[i] += mods1[i]) >= _base[i])
962 1
					mods[i] -= _base[i];
963
			// к следующей попытке
964 1
			if (trials != SIZE_MAX && trials-- == 0)
965 0
				return FALSE;
966
		}
967
	}
968 0
	return FALSE;
969
}
970

971 1
size_t priExtendPrime_deep(size_t l, size_t n, size_t base_count)
972
{
973 1
	const size_t m = W_OF_B(l);
974 1
	const size_t mo = O_OF_B(l);
975 1
	const size_t qr_deep = zmCreate_deep(mo);
976 1
	ASSERT(m >= n);
977 1
	return O_OF_W(m - n + 1 + m + 1 + m + 2 * base_count) + zmCreate_keep(mo) +
978 1
		utilMax(4,
979
			zzDiv_deep(m, n),
980 1
			zzMul_deep(n, m - n + 1),
981
			qr_deep,
982
			qrPower_deep(m, m, qr_deep));
983
}

Read our documentation on viewing source code .

Loading