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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.171  
changed lines
  Added in v.172

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25