/[dtapublic]/to_be_filed/sf_code/esrgnxpj/sfnthcgi0304/subfunc_pfact_18.c
ViewVC logotype

Contents of /to_be_filed/sf_code/esrgnxpj/sfnthcgi0304/subfunc_pfact_18.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 29 - (show annotations) (download)
Sat Oct 8 07:08:47 2016 UTC (7 years, 5 months ago) by dashley
File MIME type: text/plain
File size: 17155 byte(s)
Directories relocated.
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 //********************************************************************************

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25