1 |
//$Header: /cvsroot/esrg/sfesrg/esrgnxpj/sfnthcgi0304/subfunc_pfact_18.c,v 1.7 2003/07/01 03:46:58 dtashley Exp $
|
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 |
//********************************************************************************
|