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

Contents of /projs/dtats/trunk/projs/2018/20180707_cgi_web_tools_aux_exe/dtats_cgi_aux_arith_large/subfunc_pfact_18.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 207 - (show annotations) (download)
Sun Jul 15 21:50:56 2018 UTC (6 years, 4 months ago) by dashley
File MIME type: text/plain
File size: 16545 byte(s)
Change keyword/EOL properties.
Correct tab characters and end-of-line whitespace.
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 //********************************************************************************

Properties

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

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25