1 |
//$Header$ |
2 |
// |
3 |
//******************************************************************************** |
4 |
//Copyright (C) 2003 David T. Ashley |
5 |
//******************************************************************************** |
6 |
//This program or source file is free software; you can redistribute it and/or |
7 |
//modify it under the terms of the GNU General Public License as published by |
8 |
//the Free Software Foundation; either version 2 of the License, or (at your |
9 |
//option) any later version. |
10 |
// |
11 |
//This program or source file is distributed in the hope that it will |
12 |
//be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 |
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 |
//GNU General Public License for more details. |
15 |
// |
16 |
//You may have received a copy of the GNU General Public License |
17 |
//along with this program; if not, write to the Free Software |
18 |
//Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
19 |
//******************************************************************************** |
20 |
// |
21 |
//This module attempts to factor an integer of up to 18 decimal |
22 |
//digits using a sieve method and subject to elapsed time constraints. |
23 |
//The constraint of "up to 18 digits" |
24 |
//is because that way trial divisors can stay in one unsigned |
25 |
//integer, which speeds the division. Additionally, there isn't |
26 |
//much use in trying to factor larger integers on a web page. |
27 |
//This subfunction is "parlor trivia" grade only. |
28 |
// |
29 |
//INPUT PARAMETERS |
30 |
//---------------- |
31 |
//This subfunction accepts the following parameters, in order. |
32 |
// |
33 |
// (a) The number to factor. |
34 |
// |
35 |
// (b) The number of Miller-Rabin rounds to use when looking |
36 |
// at the probable primality of the original argument or of |
37 |
// residues. |
38 |
// |
39 |
// (c) The maximum time that should be allowed to elapse before the |
40 |
// program terminates (perhaps without finding factors). |
41 |
// This constraint is to protect the server capacity. (In the long |
42 |
// term, however, the web page may need to be removed--if too many |
43 |
// people use it, it will bring any server to a crawl.) |
44 |
// |
45 |
//OUTPUT RESULTS |
46 |
//-------------- |
47 |
//The notation below gives the outputs of the program. In some cases, [i] notation |
48 |
//is used to indicate line numbers. |
49 |
// |
50 |
//[01] An overall success or failure code for the operation. Valid responses |
51 |
// are: |
52 |
// S : Success (for this subfunction, the only possible outcome). |
53 |
// |
54 |
//[02] The fully normalized first integer entered. This means it has been |
55 |
// stripped of all weird characters, etc., and also perhaps assigned |
56 |
// a default value if the original parameter wasn't acceptable. |
57 |
// This can be used by the PHP script to repopulate form boxes. |
58 |
//[03] The fully normalized second integer entered. |
59 |
//[04] The fully normalized second integer entered. |
60 |
//[05] The type of factor specified on the next line. |
61 |
// Possibilities are: |
62 |
// "P" : The factor is definitely prime. |
63 |
// "C" : The factor is definitely composite. |
64 |
// "p" : The factor is probably prime, established by |
65 |
// Miller-Rabin. |
66 |
// |
67 |
// Note that the "C" code can only appear if the utility has to give |
68 |
// up because it runs out of allowed time. This means that the input |
69 |
// number has not been fully factored or perhaps not factored at all. |
70 |
// The "C" code can only appear for the last factor or the only factor. |
71 |
// |
72 |
// The "p" code can also only appear for the last factor. This occurs when |
73 |
// either the original input argument is prime as established by |
74 |
// Miller-Rabin or else a division leaves a result that is similarly |
75 |
// established by Miller-Rabin. |
76 |
//[06] The factor itself. |
77 |
//[07] Its multiplicity. |
78 |
//[..] Lines 5-7 are repeated for as many factors as are located. |
79 |
//[..] The second-to-last line will contain "X". |
80 |
//[..] The last line will contain an "S". |
81 |
// |
82 |
//Note the following: |
83 |
// (a) Any valid output will contain at least 9 lines. |
84 |
// (b) Any valid output will have a number of lines divisible by 3. |
85 |
// (c) The number of factors found is (nlines - 6)/3 or alternately |
86 |
// nlines/3 - 2. |
87 |
// (d) What happened can be determined using the number of lines plus |
88 |
// the code of the last factor. |
89 |
// |
90 |
//The return value (exit code) from this subfunction is always 0. |
91 |
// |
92 |
|
93 |
#define MODULE_SUBFUNC_PFACT_18 |
94 |
|
95 |
#include <assert.h> |
96 |
#include <ctype.h> |
97 |
#include <stddef.h> |
98 |
#include <stdio.h> |
99 |
#include <stdlib.h> |
100 |
#include <string.h> |
101 |
#include <time.h> |
102 |
|
103 |
#include <gmp.h> |
104 |
|
105 |
#include "auxfuncs.h" |
106 |
#include "subfunc_pfact_18.h" |
107 |
#include "sieve_eratosthenes.h" |
108 |
|
109 |
|
110 |
int SUBFUNC_PFACT_18_main(int argc, char *argv[]) |
111 |
{ |
112 |
//The time snapshot against which we compare to see if we're over |
113 |
//time budget. |
114 |
time_t time_snapshot; |
115 |
|
116 |
//Normalized first and second parameters (the integers). |
117 |
char *arg1; |
118 |
char *arg2; |
119 |
char *arg3; |
120 |
|
121 |
//Number of Miller-Rabin iterations to establish probable primality. |
122 |
int miller_rabin_iterations; |
123 |
|
124 |
//Maximum elapsed time allowed. |
125 |
int max_time; |
126 |
|
127 |
//Temporary iteration integers. |
128 |
int i; |
129 |
|
130 |
//Mask counter. Only check for timeout periodically, as calling the OS to get time |
131 |
//is presumed expensive. |
132 |
int mask_counter; |
133 |
|
134 |
//The current trial divisor. |
135 |
unsigned long trial_divisor, new_trial_divisor; |
136 |
|
137 |
//The current sieve table index. |
138 |
int sieve_table_index; |
139 |
|
140 |
//An exit flag kept to remember if we should bail the loop. |
141 |
int exit_flag; |
142 |
|
143 |
//The number to factor. |
144 |
mpz_t number_to_factor; |
145 |
|
146 |
//The square root limit. We only need to go that far. |
147 |
mpz_t square_root_limit; |
148 |
|
149 |
//The multiplicity of any factors we find. |
150 |
int multiplicity; |
151 |
|
152 |
//Initialize all of the GMP variables. |
153 |
mpz_init(number_to_factor); |
154 |
mpz_init(square_root_limit); |
155 |
|
156 |
//Grab a time snapshot. |
157 |
time_snapshot = time(NULL); |
158 |
|
159 |
//There must be an acceptable number of command-line arguments. If not, |
160 |
//abort the progam with phony data. |
161 |
if (argc != 5) |
162 |
{ |
163 |
printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); |
164 |
return(0); |
165 |
} |
166 |
|
167 |
//Copy the command-line arguments to a safe place where we can manipulate them. |
168 |
//Leave 2 characters of space in case we assign a "0". |
169 |
arg1 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[2])) + 1) * sizeof(char)); |
170 |
arg2 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[3])) + 1) * sizeof(char)); |
171 |
arg3 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[4])) + 1) * sizeof(char)); |
172 |
if ((arg1 == NULL) || (arg2 == NULL) || (arg3 == NULL)) |
173 |
{ |
174 |
printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); |
175 |
return(0); |
176 |
} |
177 |
strcpy(arg1, argv[2]); |
178 |
strcpy(arg2, argv[3]); |
179 |
strcpy(arg3, argv[4]); |
180 |
|
181 |
//Strip all of the non-digit trash out of the arguments. |
182 |
AUXFUNCS_remove_non_digits(arg1); |
183 |
AUXFUNCS_remove_non_digits(arg2); |
184 |
AUXFUNCS_remove_non_digits(arg3); |
185 |
|
186 |
//Remove all leading zeros from arguments. |
187 |
AUXFUNCS_remove_leading_zeros(arg1); |
188 |
AUXFUNCS_remove_leading_zeros(arg2); |
189 |
AUXFUNCS_remove_leading_zeros(arg3); |
190 |
|
191 |
//If an argument is zero length, fill it in with 0. |
192 |
if (!strlen(arg1)) |
193 |
strcpy(arg1, "0"); |
194 |
if (!strlen(arg2)) |
195 |
strcpy(arg2, "0"); |
196 |
if (!strlen(arg3)) |
197 |
strcpy(arg3, "0"); |
198 |
|
199 |
//We are not allowed to have 0's in this function, so |
200 |
//abort on zeros. Also, we can't have 1 for a number to check, |
201 |
//as 1 can't be factored. |
202 |
if ((!strcmp(arg1, "0")) || (!strcmp(arg2, "0")) || (!strcmp(arg3, "0")) || (!strcmp(arg1, "1"))) |
203 |
{ |
204 |
printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); |
205 |
return(0); |
206 |
} |
207 |
|
208 |
//If the number to factor exceeds 18 digits, abort. Anything that has |
209 |
//18 or fewer digits is allowed. |
210 |
if (strlen(arg1) >18) |
211 |
{ |
212 |
printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); |
213 |
return(0); |
214 |
} |
215 |
|
216 |
//Assign the number to factor. This is definitely a |
217 |
//GMP type. |
218 |
mpz_set_str(number_to_factor, arg1, 10); |
219 |
|
220 |
//Assign the number of Miller-Rabin repetitons to use, subject to an |
221 |
//absolute maximum of 1000 and minimum of 1. |
222 |
miller_rabin_iterations = 25; |
223 |
if (strlen(arg2) > 3) |
224 |
{ |
225 |
miller_rabin_iterations = 1000; |
226 |
} |
227 |
else |
228 |
{ |
229 |
sscanf(arg2, "%d", &miller_rabin_iterations); |
230 |
if (miller_rabin_iterations < 1) |
231 |
miller_rabin_iterations = 1; |
232 |
} |
233 |
|
234 |
//Assign the maximum time allowed, subject to a maximum of 1000 and minimum of 3. |
235 |
if (strlen(arg3) > 3) |
236 |
{ |
237 |
max_time = 1000; |
238 |
} |
239 |
else |
240 |
{ |
241 |
sscanf(arg3, "%d", &max_time); |
242 |
if (max_time < 3) |
243 |
max_time = 3; |
244 |
} |
245 |
|
246 |
//Output the header information before beginning the search. |
247 |
printf("S\n"); |
248 |
mpz_out_str(stdout, 10, number_to_factor); |
249 |
printf("\n"); |
250 |
printf("%d\n", miller_rabin_iterations); |
251 |
printf("%d\n", max_time); |
252 |
|
253 |
//Loop through the list of divisors we should try before starting the sieve. |
254 |
//This won't possibly exceed our time, so do it. |
255 |
for (i=0; i<SIEVE_ERATOSTHENES_N_SIEVE_FACTORS; i++) |
256 |
{ |
257 |
multiplicity = 0; |
258 |
|
259 |
//printf("Trial divisor: %d\n", SIEVE_ERATOSTHENES_sieve_factors[i]); |
260 |
|
261 |
//Factor out all occurrences that we can. |
262 |
while (mpz_divisible_ui_p(number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i])) |
263 |
{ |
264 |
mpz_divexact_ui(number_to_factor, number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i]); |
265 |
multiplicity++; |
266 |
} |
267 |
|
268 |
//Now output the information record to the stdout if we could factor it. |
269 |
if (multiplicity) |
270 |
{ |
271 |
printf("P\n"); |
272 |
printf("%u\n", SIEVE_ERATOSTHENES_sieve_factors[i]); |
273 |
printf("%d\n", multiplicity); |
274 |
} |
275 |
} |
276 |
|
277 |
//Gear up for tabulated operation as we sieve. |
278 |
new_trial_divisor = trial_divisor = SIEVE_ERATOSTHENES_FIRST_TRIAL_DIVISOR; |
279 |
sieve_table_index = SIEVE_ERATOSTHENES_FIRST_SIEVE_INDEX; |
280 |
mask_counter = 0; |
281 |
|
282 |
//Set the exit flag initially. One thing that could have |
283 |
//caused us to be ready to exit already is if we brought the |
284 |
//number to factor down to 1. We check this now because |
285 |
//otherwise only check it when have factored something out and |
286 |
//seive loop would run indefinitely if didn't do this. |
287 |
exit_flag = 0; |
288 |
if (!mpz_cmp_ui(number_to_factor, 1)) |
289 |
exit_flag = 1; |
290 |
|
291 |
//Effectively, we may have broken down the number to factor |
292 |
//with just a few operations. So, effectively, we are |
293 |
//starting afresh here. We only need to proceed to the |
294 |
//square root of the current number to factor (not the |
295 |
//original one. |
296 |
mpz_sqrt(square_root_limit, number_to_factor); |
297 |
|
298 |
//If the last value from the numbers to try before we start |
299 |
//the sieve is already past the square root limit, then whatever |
300 |
//remains is definitely prime, and we can exit. |
301 |
if (mpz_cmp_ui(square_root_limit, SIEVE_ERATOSTHENES_sieve_factors[SIEVE_ERATOSTHENES_N_SIEVE_FACTORS - 1]) <= 0) |
302 |
{ |
303 |
if (mpz_cmp_ui(number_to_factor, 1)) |
304 |
{ |
305 |
exit_flag = 1; |
306 |
printf("P\n"); |
307 |
mpz_out_str(stdout, 10, number_to_factor); |
308 |
printf("\n"); |
309 |
printf("1\n"); |
310 |
goto done; |
311 |
} |
312 |
} |
313 |
|
314 |
//If Miller-Rabin says that the remaining number is probably prime, that is good |
315 |
//enough. |
316 |
if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1) |
317 |
{ |
318 |
exit_flag = 1; |
319 |
printf("p\n"); |
320 |
mpz_out_str(stdout, 10, number_to_factor); |
321 |
printf("\n"); |
322 |
printf("1\n"); |
323 |
goto done; |
324 |
} |
325 |
|
326 |
//Loop entry condition. |
327 |
multiplicity = 0; |
328 |
|
329 |
//This is the main loop. We only check for elapsed time |
330 |
//infrequently. |
331 |
while (!exit_flag) |
332 |
{ |
333 |
//Replace the trial divisor. |
334 |
trial_divisor = new_trial_divisor; |
335 |
|
336 |
//printf("Trial divisor: %d\n", trial_divisor); |
337 |
|
338 |
//Factor out all occurrences that we can of the trial divisor, if we can. |
339 |
while (mpz_divisible_ui_p(number_to_factor, trial_divisor)) |
340 |
{ |
341 |
mpz_divexact_ui(number_to_factor, number_to_factor, trial_divisor); |
342 |
multiplicity++; |
343 |
} |
344 |
|
345 |
//Now output the information record to the stdout if we could factor it. |
346 |
if (multiplicity) |
347 |
{ |
348 |
printf("P\n"); |
349 |
printf("%u\n", trial_divisor); |
350 |
printf("%d\n", multiplicity); |
351 |
|
352 |
//We need to calculate a new square root bound. |
353 |
mpz_sqrt(square_root_limit, number_to_factor); |
354 |
|
355 |
//If we are over the square root bound, the remaining number to factor |
356 |
//is prime and we should exit. |
357 |
if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) |
358 |
exit_flag = 1; |
359 |
|
360 |
//Are we over the time bound? We might has well check here, |
361 |
//because finding prime factors is rare. |
362 |
if ((time(NULL) - time_snapshot) > max_time) |
363 |
exit_flag = 1; |
364 |
|
365 |
//Does the Miller-Rabin test say that the remaining number |
366 |
//is with near perfect certainty prime? |
367 |
if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) >= 1) |
368 |
exit_flag = 1; |
369 |
|
370 |
//Is the remaining number "1", indicating we've factored it fully? |
371 |
if (!mpz_cmp_ui(number_to_factor, 1)) |
372 |
exit_flag = 1; |
373 |
|
374 |
//Set multiplicity to zero again so don't re-enter until successful |
375 |
//division again. |
376 |
multiplicity = 0; |
377 |
} |
378 |
|
379 |
//Advance to our next trial divisor. |
380 |
new_trial_divisor = trial_divisor + SIEVE_ERATOSTHENES_sieve[sieve_table_index]; |
381 |
|
382 |
//If the new is < the old, we've rolled over. This means an exit is necessary. |
383 |
if (new_trial_divisor < trial_divisor) |
384 |
exit_flag = 1; |
385 |
|
386 |
//Advance the sieve index. |
387 |
sieve_table_index = (sieve_table_index + 1) % SIEVE_ERATOSTHENES_N_SIEVE; |
388 |
|
389 |
//This is our only chance to check for termination conditions that |
390 |
//don't come about from a successful division. But we don't do |
391 |
//this often. |
392 |
mask_counter++; |
393 |
if (!(mask_counter & 0xFFFFF)) |
394 |
{ |
395 |
//First, have we exceeded the square root bound? |
396 |
if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) |
397 |
exit_flag = 1; |
398 |
//Second, are we over the time budget? |
399 |
if ((time(NULL) - time_snapshot) > max_time) |
400 |
exit_flag = 1; |
401 |
} |
402 |
} |
403 |
|
404 |
//If we've made it out of the loop, there are a variety of reasons for that. |
405 |
//Find the right one and close up. |
406 |
if (!mpz_cmp_ui(number_to_factor, 1)) |
407 |
{ |
408 |
//We divided the number down to 1. There is nothing further to do. |
409 |
} |
410 |
else if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) |
411 |
{ |
412 |
//We are at or over the square root limit. The remaining number is definitely prime. |
413 |
printf("P\n"); |
414 |
mpz_out_str(stdout, 10, number_to_factor); |
415 |
printf("\n"); |
416 |
printf("1\n"); |
417 |
} |
418 |
else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 0) |
419 |
{ |
420 |
//Miller-Rabin says the remaining number is definitely composite. |
421 |
printf("C\n"); |
422 |
mpz_out_str(stdout, 10, number_to_factor); |
423 |
printf("\n"); |
424 |
printf("1\n"); |
425 |
} |
426 |
else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1) |
427 |
{ |
428 |
//Miller-Rabin says the remaining number is probably prime. |
429 |
printf("p\n"); |
430 |
mpz_out_str(stdout, 10, number_to_factor); |
431 |
printf("\n"); |
432 |
printf("1\n"); |
433 |
} |
434 |
else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 2) |
435 |
{ |
436 |
//Miller-Rabin says the remaining number is definitely prime. |
437 |
printf("P\n"); |
438 |
mpz_out_str(stdout, 10, number_to_factor); |
439 |
printf("\n"); |
440 |
printf("1\n"); |
441 |
} |
442 |
|
443 |
//Output the invariant footer information. |
444 |
done: |
445 |
printf("X\n"); |
446 |
printf("S\n"); |
447 |
|
448 |
//Always return 0. |
449 |
return(0); |
450 |
} |
451 |
|
452 |
//******************************************************************************** |
453 |
// $Log: subfunc_pfact_18.c,v $ |
454 |
// Revision 1.7 2003/07/01 03:46:58 dtashley |
455 |
// Edits towards working continued fraction best rational approximation |
456 |
// functionality. |
457 |
// |
458 |
// Revision 1.6 2003/04/17 20:02:05 dtashley |
459 |
// License text for the GPL added. All source files are now under the GPL, |
460 |
// after some discussion on the GMP list. |
461 |
// |
462 |
// Revision 1.5 2003/04/16 07:22:37 dtashley |
463 |
// All checks completed. |
464 |
// |
465 |
// Revision 1.4 2003/04/16 07:02:06 dtashley |
466 |
// Spelling of Greek name corrected to Eratosthenes from incorrect Erastothenes. |
467 |
// |
468 |
// Revision 1.3 2003/04/16 06:49:21 dtashley |
469 |
// Edits. |
470 |
// |
471 |
// Revision 1.2 2003/04/16 03:25:15 dtashley |
472 |
// Seems to be working correctly. Only a careful proofreading and some |
473 |
// testing remain. |
474 |
// |
475 |
// Revision 1.1 2003/04/15 23:54:58 dtashley |
476 |
// Initial checkin. |
477 |
//******************************************************************************** |
478 |
// End of SUBFUNC_PFACT_18.C. |
479 |
//******************************************************************************** |