/[dtapublic]/projs/dtats/trunk/projs/2018/20180707_cgi_web_tools_aux_exe/subfunc_pfact_18.c
ViewVC logotype

Annotation of /projs/dtats/trunk/projs/2018/20180707_cgi_web_tools_aux_exe/subfunc_pfact_18.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 172 - (hide annotations) (download)
Sun Jul 8 21:17:26 2018 UTC (5 years, 7 months ago) by dashley
File MIME type: text/plain
File size: 16573 byte(s)
Add keyword expansion, native EOL.
1 dashley 172 //$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     //********************************************************************************

Properties

Name Value
svn:eol-style native
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25