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

Properties

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

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25