2018/20180707_cgi_web_tools_aux_exe/subfunc_pfact_18.c

//$Header: /cvsroot/esrg/sfesrg/esrgnxpj/sfnthcgi0304/subfunc_pfact_18.c,v 1.7 2003/07/01 03:46:58 dtashley Exp $
//
//********************************************************************************
//Copyright (C) 2003 David T. Ashley
//********************************************************************************
//This program or source file is free software; you can redistribute it and/or 
//modify it under the terms of the GNU General Public License as published by
//the Free Software Foundation; either version 2 of the License, or (at your 
//option) any later version.
//
//This program or source file is distributed in the hope that it will 
//be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.
//
//You may have received a copy of the GNU General Public License
//along with this program; if not, write to the Free Software
//Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//********************************************************************************
//
//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);
   } 

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