prime.d

Revision 28, 12.8 kB (checked in by Arcane Jill, 5 years ago)
--

Line
1	module etc.bigint.prime;
2	import etc.bigint.bigint_int;
3	import etc.bigint.modexp;
4	import etc.prime;
5	private import std.random;
6
7	/*
8
9	Copyright (c) 2004, Arcane Jill
10
11	All rights reserved. Intellectual Property Me Arse!
12
13	Redistribution and use in source and binary forms, with or without modification, are permitted
14	provided that the following conditions are met:
15
16	* Redistributions of source code must retain the above copyright notice, the phrase
17	"Intellectual Property Me Arse!", this list of conditions, and the following disclaimer.
18	* Redistributions in binary form must reproduce the above copyright notice, the phrase
19	"Intellectual Property Me Arse!", this list of conditions and the following disclaimer
20	in the documentation and/or other materials provided with the distribution.
21	* The name Arcane Jill may not be used to endorse or promote products derived from this
22	software without specific prior written permission.
23
24	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS
25	OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
26	AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER,
27	OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
29	SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED, AND ON ANY
30	THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
31	OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32	OF SUCH DAMAGE.
33
34	*/
35
36	/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37	This is the all-important primality test
38	*/
39
40	enum { DEFINITELY_FALSE=0, PROBABLY_TRUE, DEFINITELY_TRUE };
41
42	int isProbablyPrime(Int p)
43	{
44	return isProbablyPrime(p, 10);
45	}
46
47	int isProbablyPrime(Int p, uint attempts)
48	out(r)
49	{
50	switch (r)
51	{
52	case DEFINITELY_FALSE:
53	case PROBABLY_TRUE:
54	case DEFINITELY_TRUE:
55	break;
56	default:
57	assert(false);
58	}
59	version(Paranoid)
60	{
61	// In the Paranoid build, we confirm that, for numbers < 0x100000000,
62	// isPrime() and isProbablyPrime(), both agree with each other.
63
64	if (p.isUint())
65	{
66	uint n = p.toUint();
67	assert ((r!=0) == isPrime(n));
68	}
69	}
70	}
71	body
72	{
73	if (p < 2) return DEFINITELY_FALSE; // This takes care of 1, 0 and all negative numbers
74	// What we are left with is guaranteed positiive
75	if (p.isUint())
76	{
77	uint n = p.toUint();
78	{
79	version(Paranoid)
80	{
81	if (n < 4) return DEFINITELY_TRUE; // In a paranoid build, we remain "unsure" about any
82	// number bigger than 4.
83	}
84	else
85	{
86	// In a non-paranoid build, we can be certain about every uint.
87	return isPrime(n) ? DEFINITELY_TRUE : DEFINITELY_FALSE;
88	}
89	}
90	}
91	if ((p & 1u) == 0)
92	{
93	return DEFINITELY_FALSE;
94	}
95	for (int i=attempts-1; i>=0; --i)
96	{
97	if (!isProbablyPrimeLegendre(p))
98	{
99	return DEFINITELY_FALSE;
100	}
101	}
102	return PROBABLY_TRUE;
103	}
104
105	private
106	{
107	bool isProbablyPrimeLegendre(Int p)
108	in
109	{
110	assert(p > 4);
111	assert((p & 1u) == 1);
112	}
113	body
114	{
115	Int a;
116	uint[] z;
117	z.length = p.end;
118	uint mask = (1 << (bsr(p[p.end-1]) + 1)) - 1;
119	do
120	{
121	for (int i=0; i<p.end; ++i)
122	{
123	z[i] = rand(); // OK, because this does not need to be cryptographically strong
124	}
125	z[p.end-1] &= mask;
126	a = new Int(z, false);
127	}
128	while (a >= p);
129
130	int j = jacobi(a,p);
131	Int lhs = (j<0) ? p+j : new Int(j);
132	Int e = (p-1)>>>1;
133	Int rhs = modPow(a,e,p);
134	return (lhs == rhs) ? true : false;
135	}
136	}
137
138	/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139	Functions to find nearby primes
140	*/
141
142	Int getProbablePrimeGreaterEqual(Int n)
143	{
144	if (n <= 2) return Int.TWO;
145	for (n=n\|1; n!=1; n=n+2)
146	{
147	if (isProbablyPrime(n)) return n;
148	}
149	}
150
151	Int getProbablePrimeGreaterEqual(Int n, uint attempts)
152	{
153	if (n <= 2) return Int.TWO;
154	for (n=n\|1; n!=1; n=n+2)
155	{
156	if (isProbablyPrime(n, attempts)) return n;
157	}
158	}
159
160	Int getProbablePrimeLessEqual(Int n)
161	{
162	if (n < 2) throw new IntException("getProbablePrimeLessEqual(n) not defined for n < 2");
163	for (n=(n\|1)-2; n!=1; n=n-2)
164	{
165	if (isProbablyPrime(n)) return n;
166	}
167	return Int.TWO;
168	}
169
170	Int getProbablePrimeLessEqual(Int n, uint attempts)
171	{
172	if (n < 2) throw new IntException("getProbablePrimeLessEqual(n) not defined for n < 2");
173	for (n=(n\|1)-2; n!=1; n=n-2)
174	{
175	if (isProbablyPrime(n, attempts)) return n;
176	}
177	return Int.TWO;
178	}
179
180	Int getProbablePrimeGreater(Int n)
181	{
182	return getProbablePrimeGreaterEqual(n+1);
183	}
184
185	Int getProbablePrimeGreater(Int n, uint attempts)
186	{
187	return getProbablePrimeGreaterEqual(n+1, attempts);
188	}
189
190	Int getProbablePrimeLess(Int n)
191	{
192	return getProbablePrimeLessEqual(n-1);
193	}
194
195	Int getProbablePrimeLess(Int n, uint attempts)
196	{
197	return getProbablePrimeLessEqual(n-1, attempts);
198	}
199
200	/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201	The Jacobi symbol
202	*/
203
204	int jacobi(Int a, Int n)
205	in
206	{
207	assert(n >= 0);
208	assert((n & 1u) == 1);
209	}
210	body
211	{
212	Int b;
213	int j = 1;
214	uint n8;
215	while (!a.equalsZero)
216	{
217	while ((a & 1u) == 0)
218	{
219	a = a >>> 1u;
220	n8 = n & 7u;
221	if ((n8 == 3u) \|\| (n8 == 5u)) j = -j;
222	}
223	b = a;
224	a = n;
225	n = b;
226	if (((a & 3u) == 3u) && ((n & 3u) == 3u)) j = -j;
227	a = a % n;
228	}
229	return (n == 1) ? j : 0;
230	}
231
232	//=======================================================================================================
233
234	version(Paranoid)
235	{
236	unittest
237	{
238	// Give the jacobi function a damn good test
239	Int a,b;
240
241	a = new Int("0xB491087F_5E666C48_96D44FDE_E6EEB574");
242	b = new Int("0x74FB8A89_127E3A79_3622C8FB_85CD2847");
243	assert(jacobi(a,b) == 1);
244	a = new Int("0x4AF1671C_077E3FB3_68FF7D90_3D953866");
245	b = new Int("0x41971D2A_8E282742_C52C8E37_064B75CB");
246	assert(jacobi(a,b) == -1);
247	a = new Int("0x5B8512B9_4A5B7F70_A455910B_ECAF73E3");
248	b = new Int("0x8D0B977C_FF670D67_56E7826C_D9DC3A5D");
249	assert(jacobi(a,b) == 1);
250	a = new Int("0x4C49A1A1_F4048915_7E799556_6E1BB14B");
251	b = new Int("0xCD0763B4_0D7394D3_2474EE00_D11D632B");
252	assert(jacobi(a,b) == -1);
253	a = new Int("0x5F559504_2F18DB9D_F18E5D5E_1195F82A");
254	b = new Int("0xB2B0286A_DA9C4BE6_D3F12D6B_410D693F");
255	assert(jacobi(a,b) == 1);
256	a = new Int("0xF5037185_7377BB04_D92E30E3_B6E3C9FD");
257	b = new Int("0xFC769BDF_5CBB5AA7_3E686320_EB58403B");
258	assert(jacobi(a,b) == 1);
259	a = new Int("0xB2A1C094_785D3393_79F7BC91_10098651");
260	b = new Int("0x108E426F_9FBA3AF1_37031C6F_453DAFD7");
261	assert(jacobi(a,b) == -1);
262	a = new Int("0xE96A98B0_DD5E858A_25FA7404_C02E85AB");
263	b = new Int("0x96803860_A7337003_E2070241_47E0D119");
264	assert(jacobi(a,b) == -1);
265	a = new Int("0xA7596AFD_5F8B5115_6B9A2110_FDE88145");
266	b = new Int("0x1A66DF32_D4A14974_DDBD1217_D2A07F59");
267	assert(jacobi(a,b) == 0);
268	a = new Int("0xA7B9B5E7_A8B0BC33_BDE83F60_F66578D0");
269	b = new Int("0x05A4028B_83740BEF_BDCD487D_73BD7169");
270	assert(jacobi(a,b) == 1);
271	a = new Int("0xF6C77351_9C893E0D_C21B778A_15E098BA");
272	b = new Int("0xF91A1784_DD65B896_DE50D24D_586B30F5");
273	assert(jacobi(a,b) == 0);
274	a = new Int("0x2B40C766_42736D1F_729AE1A6_CAF09121");
275	b = new Int("0x51093349_B19E34EE_26410AC2_C5B3EF31");
276	assert(jacobi(a,b) == 1);
277	a = new Int("0x01B33E51_F6430CDF_159785FF_B6EECADF");
278	b = new Int("0x6C0F696F_EEDE49F4_00D3D17C_E421C727");
279	assert(jacobi(a,b) == 1);
280	a = new Int("0x762F792F_FF57DF2B_2B49E8C5_6FCE7908");
281	b = new Int("0xA835C636_64CEB4D5_66505CB8_9E9D7603");
282	assert(jacobi(a,b) == -1);
283	a = new Int("0x4722A982_CCE680C9_18FDA4EE_EEDB4F03");
284	b = new Int("0xF99BF268_C4AE0931_73BD44FF_2A09334F");
285	assert(jacobi(a,b) == 0);
286	a = new Int("0x9DE41532_EE4ED214_36CB2F10_92F40D28");
287	b = new Int("0xAFAB4F44_DE5248A9_B33910C6_E2872B91");
288	assert(jacobi(a,b) == 1);
289	a = new Int("0x8472F4D0_A030C25D_DCA5E634_54316504");
290	b = new Int("0x2E0A2913_8709D09C_A99B993B_7411F00D");
291	assert(jacobi(a,b) == 1);
292	a = new Int("0x71E4E601_EFB6F31D_CB6EB1A0_80145940");
293	b = new Int("0x799FB97E_833EDAE2_B2FFE313_DD15C29B");
294	assert(jacobi(a,b) == 1);
295	a = new Int("0x89E95A10_EA7F090F_89643BCC_E1CBA334");
296	b = new Int("0x9EFCF2F0_E9A5EA86_38F6DDCF_70D072C3");
297	assert(jacobi(a,b) == -1);
298	a = new Int("0x3553F530_A3CCE32B_1417D3AC_FAF4AA50");
299	b = new Int("0x62B4B624_F2D6A117_E8AAA16A_4DAD2857");
300	assert(jacobi(a,b) == 1);
301	a = new Int("0x78BF9779_D01631B6_57FB1F75_BB503B05");
302	b = new Int("0xD2EE41AC_0F267F82_B608D8BD_5EAEC2A7");
303	assert(jacobi(a,b) == 1);
304	a = new Int("0xBEC7C504_662EAC88_110F8E7F_41CF6716");
305	b = new Int("0x9B2867AF_8B04E9DC_3A2CB599_FF88AC1B");
306	assert(jacobi(a,b) == -1);
307	a = new Int("0x88EBA6DC_68AE8A25_4F982831_ACD3EE37");
308	b = new Int("0xB3B110DE_CF8AE31C_A2762C5B_DA7C2213");
309	assert(jacobi(a,b) == 1);
310	a = new Int("0xBDA16A3A_21924BF8_8E4615B1_AEC224A6");
311	b = new Int("0x63071AC6_D101CC44_9C311E22_C21B1325");
312	assert(jacobi(a,b) == 1);
313	a = new Int("0xE7E37905_470BCEC5_387A2281_80722136");
314	b = new Int("0x164C6650_3F788811_837B2BDC_7C81E791");
315	assert(jacobi(a,b) == -1);
316	a = new Int("0x06058196_FDFEC957_FF49624E_8F8046B6");
317	b = new Int("0x1E3232DC_1F78C8E3_CC9B01FC_2C7E3387");
318	assert(jacobi(a,b) == 0);
319	a = new Int("0x344F20C1_95C1666C_6FE4E8E1_CF7148DD");
320	b = new Int("0xB0927117_21B633E4_6D827CC0_7B5F1921");
321	assert(jacobi(a,b) == 1);
322	a = new Int("0xAD87C698_08FB6D1E_D38E6697_5BE4DE28");
323	b = new Int("0x522CDA03_19750539_30AD4238_CB2F4771");
324	assert(jacobi(a,b) == -1);
325	a = new Int("0xC8F8D9B0_CB82E686_F6FA9581_A09E278F");
326	b = new Int("0xF37979EF_062ED557_E5B2F0D2_CCA080CF");
327	assert(jacobi(a,b) == 1);
328	a = new Int("0xFDA99FCB_74E6A156_052AD4AE_4608B8E2");
329	b = new Int("0xC1713C5A_EFD30183_E7877937_5A6C809F");
330	assert(jacobi(a,b) == 1);
331	a = new Int("0x5FCB6751_EE7929E3_FB85A93A_D74F9012");
332	b = new Int("0x21739F10_3E8DB626_2CD355FA_83334A6F");
333	assert(jacobi(a,b) == -1);
334	a = new Int("0x09D8AF93_19F3EA16_1FF1DFB4_0F498670");
335	b = new Int("0x3CFEFC1E_962526D5_F9319E3C_9B942461");
336	assert(jacobi(a,b) == -1);
337	a = new Int("0x75FD6803_868338E0_403E9FAC_38ED7AC1");
338	b = new Int("0xC4505E5A_A3D1F623_6B481B80_B71E78A1");
339	assert(jacobi(a,b) == 1);
340	a = new Int("0xC8C86E65_C0C58EEE_EC63AEE2_B2682CBA");
341	b = new Int("0x2EACE0BF_BC270229_5767F1A0_1A22A209");
342	assert(jacobi(a,b) == -1);
343	a = new Int("0xC57B8E02_EA42E821_CA01A42F_8C618E1F");
344	b = new Int("0x8573627C_F0FD8735_29E580A6_12A9E975");
345	assert(jacobi(a,b) == -1);
346	a = new Int("0x22C5BEFE_A1364B1E_7B614CF2_E4FADF75");
347	b = new Int("0x770000E8_D3ED7103_1B497D5C_D3969B81");
348	assert(jacobi(a,b) == -1);
349	a = new Int("0x4A279993_6B345BB8_7217CE93_82A33EC9");
350	b = new Int("0x1CC25ECE_641199F5_C5788C02_D61F9A31");
351	assert(jacobi(a,b) == 0);
352	a = new Int("0x4727BCBF_A9C3FD02_C27F9710_E5F1F039");
353	b = new Int("0x80BF9919_EC343F5C_2CD7E971_DA8ADF0D");
354	assert(jacobi(a,b) == 1);
355	a = new Int("0x0E521046_75B30630_1DB4A9D4_5B4DA6F2");
356	b = new Int("0x817B9B73_AF2194C0_7A305E46_4717EBF5");
357	assert(jacobi(a,b) == -1);
358	a = new Int("0x68B46BE1_7E692C63_2608CE5F_277F26FE");
359	b = new Int("0x78F4F199_274B688C_DB8F28C5_9D18C5E3");
360	assert(jacobi(a,b) == 0);
361	a = new Int("0x344A01E8_BAEF8366_24CFCD49_CA532A2F");
362	b = new Int("0xDAF1A2E3_AAD2D3C8_AA85B1C4_47283DAB");
363	assert(jacobi(a,b) == -1);
364	a = new Int("0xEB53DA77_67205493_66CB5A7E_0947A58E");
365	b = new Int("0x7ED21C8F_1BE47981_EF01BF9C_12DA0297");
366	assert(jacobi(a,b) == 1);
367	a = new Int("0x916BB67B_1363CBC1_0FD1C5B5_7A018EFA");
368	b = new Int("0xB4EEE085_71ABE95F_3130EC15_A83D19F3");
369	assert(jacobi(a,b) == -1);
370	a = new Int("0xB073A95E_4A431DEB_5AC7B208_D92E74D6");
371	b = new Int("0x241C1DB6_7CAE5C32_FD4D2ADB_49AE56C5");
372	assert(jacobi(a,b) == 1);
373	a = new Int("0x24795FBB_E2C87FF9_BB398AA4_73877D91");
374	b = new Int("0xD9839622_B5420184_12E8333D_CC09C5AD");
375	assert(jacobi(a,b) == 0);
376	a = new Int("0x668AAFCE_F7958D99_FD13B620_1970C8F3");
377	b = new Int("0x177063D0_7746E325_373A7774_49122A01");
378	assert(jacobi(a,b) == -1);
379	a = new Int("0x1E3D06DA_2E04958E_28042C2B_40CE928F");
380	b = new Int("0xD3E919C8_696DA631_D829E9E8_3EFB301D");
381	assert(jacobi(a,b) == -1);
382	a = new Int("0x9275F823_8E6A96AC_08DC303F_4E1608CC");
383	b = new Int("0xB3396244_95731110_4C9B5943_FEA76133");
384	assert(jacobi(a,b) == -1);
385	a = new Int("0x1F86E171_C66697AC_DD8AC83D_663F1D6F");
386	b = new Int("0x585CABB2_FD1692C3_EEC6C944_B62348F9");
387	assert(jacobi(a,b) == -1);
388	a = new Int("0xF8D6785F_659C6674_23E2D9DB_77162E9E");
389	b = new Int("0x814510FC_ECFC266F_B55C81FA_576B4B8F");
390	assert(jacobi(a,b) == 1);
391	}
392	}

Note: See TracBrowser for help on using the browser.

deimos

root/trunk/etc/bigint/prime.d

Download in other formats: