2018/20180707_cgi_web_tools_aux_exe/subfunc_pfact_18.c

//$Header$
//********************************************************************************
//Copyright (c) 2003, 2018 David T. Ashley.
//********************************************************************************
//This file is part of "arith_large_cgi", a program that is designed to be
//invoked by a PHP script as part of serving a web page that performs
//calculations involving large integers.  (A secondary compiled program is
//used because a compiled program can perform certain  calculation-intensive
//tasks far more efficiently than a PHP script.)  This program is provided by
//David T. Ashley (dashley@gmail.com) under the MIT License (reproduced
//immediately below).
//********************************************************************************
//Permission is hereby granted, free of charge, to any person obtaining a copy
//of this software and associated documentation files (the "Software"), to deal
//in the Software without restriction, including without limitation the rights
//to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
//copies of the Software, and to permit persons to whom the Software is
//furnished to do so, subject to the following conditions:
//
//The above copyright notice and this permission notice shall be included in all
//copies or substantial portions of the Software.
//
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
//SOFTWARE.
//********************************************************************************
//This module attempts to factor an integer of up to 18 decimal
//digits using a sieve method and subject to elapsed time constraints.
//The constraint of "up to 18 digits"
//is because that way trial divisors can stay in one unsigned
//integer, which speeds the division.  Additionally, there isn't
//much use in trying to factor larger integers on a web page.
//This subfunction is "parlor trivia" grade only.
// 
//INPUT PARAMETERS
//----------------
//This subfunction accepts the following parameters, in order.
//
// (a) The number to factor.
//
// (b) The number of Miller-Rabin rounds to use when looking
//     at the probable primality of the original argument or of
//     residues.
//
// (c) The maximum time that should be allowed to elapse before the
//     program terminates (perhaps without finding factors).
//     This constraint is to protect the server capacity.  (In the long
//     term, however, the web page may need to be removed--if too many
//     people use it, it will bring any server to a crawl.)
//
//OUTPUT RESULTS
//--------------
//The notation below gives the outputs of the program.  In some cases, [i] notation
//is used to indicate line numbers.
//
//[01] An overall success or failure code for the operation.  Valid responses
//     are:
//        S      :  Success (for this subfunction, the only possible outcome).
//
//[02] The fully normalized first integer entered.  This means it has been
//     stripped of all weird characters, etc., and also perhaps assigned
//     a default value if the original parameter wasn't acceptable.
//     This can be used by the PHP script to repopulate form boxes.
//[03] The fully normalized second integer entered.
//[04] The fully normalized second integer entered.
//[05] The type of factor specified on the next line.
//     Possibilities are:
//         "P"  :  The factor is definitely prime.
//         "C"  :  The factor is definitely composite.
//         "p"  :  The factor is probably prime, established by
//                 Miller-Rabin.
//
//     Note that the "C" code can only appear if the utility has to give
//     up because it runs out of allowed time.  This means that the input
//     number has not been fully factored or perhaps not factored at all.
//     The "C" code can only appear for the last factor or the only factor.
//
//     The "p" code can also only appear for the last factor.  This occurs when
//     either the original input argument is prime as established by
//     Miller-Rabin or else a division leaves a result that is similarly
//     established by Miller-Rabin.
//[06] The factor itself.
//[07] Its multiplicity.
//[..] Lines 5-7 are repeated for as many factors as are located.
//[..] The second-to-last line will contain "X".
//[..] The last line will contain an "S".
//
//Note the following:
//   (a) Any valid output will contain at least 9 lines.
//   (b) Any valid output will have a number of lines divisible by 3.
//   (c) The number of factors found is (nlines - 6)/3 or alternately
//       nlines/3 - 2.
//   (d) What happened can be determined using the number of lines plus
//       the code of the last factor.
//
//The return value (exit code) from this subfunction is always 0.
// 

#define MODULE_SUBFUNC_PFACT_18

#include <assert.h>
#include <ctype.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include <gmp.h>

#include "auxfuncs.h"
#include "subfunc_pfact_18.h"
#include "sieve_eratosthenes.h"


int SUBFUNC_PFACT_18_main(int argc, char *argv[])
   {
   //The time snapshot against which we compare to see if we're over
   //time budget.
   time_t time_snapshot;

   //Normalized first and second parameters (the integers).
   char *arg1;
   char *arg2;
   char *arg3;

   //Number of Miller-Rabin iterations to establish probable primality.
   int miller_rabin_iterations;

   //Maximum elapsed time allowed.
   int max_time;

   //Temporary iteration integers.
   int i;

   //Mask counter.  Only check for timeout periodically, as calling the OS to get time
   //is presumed expensive.
   int mask_counter;

   //The current trial divisor.
   unsigned long trial_divisor, new_trial_divisor;

   //The current sieve table index.
   int sieve_table_index;

   //An exit flag kept to remember if we should bail the loop.
   int exit_flag;

   //The number to factor.
   mpz_t number_to_factor;

   //The square root limit.  We only need to go that far.
   mpz_t square_root_limit;

   //The multiplicity of any factors we find.
   int multiplicity;

   //Initialize all of the GMP variables.
   mpz_init(number_to_factor);
   mpz_init(square_root_limit);

   //Grab a time snapshot.
   time_snapshot = time(NULL);

   //There must be an acceptable number of command-line arguments.   If not,
   //abort the progam with phony data.
   if (argc != 5)
      {
      printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n");
      return(0);
      }

   //Copy the command-line arguments to a safe place where we can manipulate them.  
   //Leave 2 characters of space in case we assign a "0".
   arg1 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[2])) + 1) * sizeof(char));
   arg2 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[3])) + 1) * sizeof(char));
   arg3 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[4])) + 1) * sizeof(char));
   if ((arg1 == NULL) || (arg2 == NULL) || (arg3 == NULL))
      {
      printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n");
      return(0);
      }
   strcpy(arg1, argv[2]);
   strcpy(arg2, argv[3]);
   strcpy(arg3, argv[4]);

   //Strip all of the non-digit trash out of the arguments.
   AUXFUNCS_remove_non_digits(arg1);
   AUXFUNCS_remove_non_digits(arg2);
   AUXFUNCS_remove_non_digits(arg3);

   //Remove all leading zeros from arguments.
   AUXFUNCS_remove_leading_zeros(arg1);
   AUXFUNCS_remove_leading_zeros(arg2);
   AUXFUNCS_remove_leading_zeros(arg3);

   //If an argument is zero length, fill it in with 0.
   if (!strlen(arg1))
      strcpy(arg1, "0");
   if (!strlen(arg2))
      strcpy(arg2, "0");
   if (!strlen(arg3))
      strcpy(arg3, "0");

   //We are not allowed to have 0's in this function, so
   //abort on zeros.  Also, we can't have 1 for a number to check,
   //as 1 can't be factored.
   if ((!strcmp(arg1, "0")) || (!strcmp(arg2, "0")) || (!strcmp(arg3, "0")) || (!strcmp(arg1, "1")))
      {
      printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n");
      return(0);
      }

   //If the number to factor exceeds 18 digits, abort.  Anything that has 
   //18 or fewer digits is allowed.
   if (strlen(arg1) >18)
      {
      printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n");
      return(0);
      }

   //Assign the number to factor.  This is definitely a
   //GMP type.
   mpz_set_str(number_to_factor, arg1, 10);

   //Assign the number of Miller-Rabin repetitons to use, subject to an 
   //absolute maximum of 1000 and minimum of 1.
   miller_rabin_iterations = 25;
   if (strlen(arg2) > 3)
      {
      miller_rabin_iterations = 1000;
      }
   else
      {
      sscanf(arg2, "%d", &miller_rabin_iterations);
      if (miller_rabin_iterations < 1)
         miller_rabin_iterations = 1;
      }

   //Assign the maximum time allowed, subject to a maximum of 1000 and minimum of 3.
   if (strlen(arg3) > 3)
      {
      max_time = 1000;
      }
   else
      {
      sscanf(arg3, "%d", &max_time);
      if (max_time < 3)
         max_time = 3;
      }

   //Output the header information before beginning the search.
   printf("S\n");
   mpz_out_str(stdout, 10, number_to_factor); 
   printf("\n");     
   printf("%d\n", miller_rabin_iterations);
   printf("%d\n", max_time);

   //Loop through the list of divisors we should try before starting the sieve.
   //This won't possibly exceed our time, so do it.
   for (i=0; i<SIEVE_ERATOSTHENES_N_SIEVE_FACTORS; i++)
      {
      multiplicity = 0;

      //printf("Trial divisor: %d\n", SIEVE_ERATOSTHENES_sieve_factors[i]);

      //Factor out all occurrences that we can.
      while (mpz_divisible_ui_p(number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i]))
         {
         mpz_divexact_ui(number_to_factor, number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i]);
         multiplicity++;
         }

      //Now output the information record to the stdout if we could factor it.
      if (multiplicity)
         {
         printf("P\n");
         printf("%u\n", SIEVE_ERATOSTHENES_sieve_factors[i]);
         printf("%d\n", multiplicity);
         }
      }
 
   //Gear up for tabulated operation as we sieve.
   new_trial_divisor = trial_divisor     = SIEVE_ERATOSTHENES_FIRST_TRIAL_DIVISOR;
   sieve_table_index                     = SIEVE_ERATOSTHENES_FIRST_SIEVE_INDEX;
   mask_counter                          = 0;

   //Set the exit flag initially.  One thing that could have
   //caused us to be ready to exit already is if we brought the
   //number to factor down to 1.  We check this now because
   //otherwise only check it when have factored something out and
   //seive loop would run indefinitely if didn't do this.
   exit_flag = 0;
   if (!mpz_cmp_ui(number_to_factor, 1))
      exit_flag = 1;

   //Effectively, we may have broken down the number to factor
   //with just a few operations.  So, effectively, we are 
   //starting afresh here.  We only need to proceed to the 
   //square root of the current number to factor (not the
   //original one.
   mpz_sqrt(square_root_limit, number_to_factor);

   //If the last value from the numbers to try before we start
   //the sieve is already past the square root limit, then whatever
   //remains is definitely prime, and we can exit.
   if (mpz_cmp_ui(square_root_limit,  SIEVE_ERATOSTHENES_sieve_factors[SIEVE_ERATOSTHENES_N_SIEVE_FACTORS - 1]) <= 0)
      {
      if (mpz_cmp_ui(number_to_factor, 1))
         {
         exit_flag = 1;
         printf("P\n");
         mpz_out_str(stdout, 10, number_to_factor); 
         printf("\n");
         printf("1\n");     
         goto done;
         }
     }

   //If Miller-Rabin says that the remaining number is probably prime, that is good
   //enough.
   if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1)
      {
      exit_flag = 1;
      printf("p\n");
      mpz_out_str(stdout, 10, number_to_factor); 
      printf("\n");
      printf("1\n");
      goto done;     
      }

   //Loop entry condition.
   multiplicity = 0;

   //This is the main loop.  We only check for elapsed time
   //infrequently.
   while (!exit_flag)
      {
      //Replace the trial divisor.
      trial_divisor = new_trial_divisor;

      //printf("Trial divisor: %d\n", trial_divisor);

      //Factor out all occurrences that we can of the trial divisor, if we can.
      while (mpz_divisible_ui_p(number_to_factor, trial_divisor))
         {
         mpz_divexact_ui(number_to_factor, number_to_factor, trial_divisor);
         multiplicity++;
         }

      //Now output the information record to the stdout if we could factor it.
      if (multiplicity)
         {
         printf("P\n");
         printf("%u\n", trial_divisor);
         printf("%d\n", multiplicity);

         //We need to calculate a new square root bound.
         mpz_sqrt(square_root_limit, number_to_factor);

         //If we are over the square root bound, the remaining number to factor
         //is prime and we should exit.
         if (mpz_cmp_ui(square_root_limit,  trial_divisor) <= 0)
            exit_flag = 1;

         //Are we over the time bound?  We might has well check here,
         //because finding prime factors is rare.
         if ((time(NULL) - time_snapshot) > max_time)
            exit_flag = 1;

         //Does the Miller-Rabin test say that the remaining number
         //is with near perfect certainty prime?
         if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) >= 1)
            exit_flag = 1;

         //Is the remaining number "1", indicating we've factored it fully?
         if (!mpz_cmp_ui(number_to_factor, 1))
            exit_flag = 1;

         //Set multiplicity to zero again so don't re-enter until successful
         //division again.
         multiplicity = 0;
         }

      //Advance to our next trial divisor.
      new_trial_divisor = trial_divisor + SIEVE_ERATOSTHENES_sieve[sieve_table_index];

      //If the new is < the old, we've rolled over.  This means an exit is necessary.
      if (new_trial_divisor < trial_divisor)
         exit_flag = 1;

      //Advance the sieve index.
      sieve_table_index = (sieve_table_index + 1) % SIEVE_ERATOSTHENES_N_SIEVE;

      //This is our only chance to check for termination conditions that
      //don't come about from a successful division.  But we don't do 
      //this often.
      mask_counter++;
      if (!(mask_counter & 0xFFFFF))
         {
         //First, have we exceeded the square root bound?
         if (mpz_cmp_ui(square_root_limit,  trial_divisor) <= 0)
            exit_flag = 1;
         //Second, are we over the time budget?
         if ((time(NULL) - time_snapshot) > max_time)
            exit_flag = 1;
         }
      }
   
   //If we've made it out of the loop, there are a variety of reasons for that.
   //Find the right one and close up.
   if (!mpz_cmp_ui(number_to_factor, 1))
      {
      //We divided the number down to 1.  There is nothing further to do.
      }
   else if (mpz_cmp_ui(square_root_limit,  trial_divisor) <= 0)
      {
      //We are at or over the square root limit.  The remaining number is definitely prime.
      printf("P\n");
      mpz_out_str(stdout, 10, number_to_factor); 
      printf("\n");
      printf("1\n");     
      }
   else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 0)
      {
      //Miller-Rabin says the remaining number is definitely composite.
      printf("C\n");
      mpz_out_str(stdout, 10, number_to_factor); 
      printf("\n");
      printf("1\n");     
      }
   else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1)
      {
      //Miller-Rabin says the remaining number is probably prime.
      printf("p\n");
      mpz_out_str(stdout, 10, number_to_factor); 
      printf("\n");
      printf("1\n");     
      }
   else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 2)
      {
      //Miller-Rabin says the remaining number is definitely prime.
      printf("P\n");
      mpz_out_str(stdout, 10, number_to_factor); 
      printf("\n");
      printf("1\n");     
      }

   //Output the invariant footer information.
   done:
   printf("X\n");
   printf("S\n");   

   //Always return 0.
   return(0);
   } 

//********************************************************************************
// End of SUBFUNC_PFACT_18.C.
//********************************************************************************
1	//$Header$
2	//********************************************************************************
3	//Copyright (c) 2003, 2018 David T. Ashley.
4	//********************************************************************************
5	//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	//
20	//The above copyright notice and this permission notice shall be included in all
21	//copies or substantial portions of the Software.
22	//
23	//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	//********************************************************************************
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	//********************************************************************************