Parent Directory | Revision Log

Revision **175** -
(**hide annotations**)
(**download**)

*Mon Jul 9 00:00:34 2018 UTC*
(6 years ago)
by *dashley*

File MIME type: text/plain

File size: 16554 byte(s)

File MIME type: text/plain

File size: 16554 byte(s)

Update license (GPL->MIT) and copyright information. Remove CVS log, which SVN does not support.

1 | dashley | 172 | //$Header$ |

2 | //******************************************************************************** | ||

3 | dashley | 175 | //Copyright (c) 2003, 2018 David T. Ashley. |

4 | dashley | 172 | //******************************************************************************** |

5 | dashley | 175 | //This file is part of "arith_large_cgi", a program that is designed to be |

6 | //invoked by a PHP script as part of serving a web page that performs | ||

7 | //calculations involving large integers. (A secondary compiled program is | ||

8 | //used because a compiled program can perform certain calculation-intensive | ||

9 | //tasks far more efficiently than a PHP script.) This program is provided by | ||

10 | //David T. Ashley (dashley@gmail.com) under the MIT License (reproduced | ||

11 | //immediately below). | ||

12 | //******************************************************************************** | ||

13 | //Permission is hereby granted, free of charge, to any person obtaining a copy | ||

14 | //of this software and associated documentation files (the "Software"), to deal | ||

15 | //in the Software without restriction, including without limitation the rights | ||

16 | //to use, copy, modify, merge, publish, distribute, sublicense, and/or sell | ||

17 | //copies of the Software, and to permit persons to whom the Software is | ||

18 | //furnished to do so, subject to the following conditions: | ||

19 | dashley | 172 | // |

20 | dashley | 175 | //The above copyright notice and this permission notice shall be included in all |

21 | //copies or substantial portions of the Software. | ||

22 | dashley | 172 | // |

23 | dashley | 175 | //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |

24 | //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | ||

25 | //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | ||

26 | //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | ||

27 | //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | ||

28 | //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE | ||

29 | //SOFTWARE. | ||

30 | dashley | 172 | //******************************************************************************** |

31 | //This module attempts to factor an integer of up to 18 decimal | ||

32 | //digits using a sieve method and subject to elapsed time constraints. | ||

33 | //The constraint of "up to 18 digits" | ||

34 | //is because that way trial divisors can stay in one unsigned | ||

35 | //integer, which speeds the division. Additionally, there isn't | ||

36 | //much use in trying to factor larger integers on a web page. | ||

37 | //This subfunction is "parlor trivia" grade only. | ||

38 | // | ||

39 | //INPUT PARAMETERS | ||

40 | //---------------- | ||

41 | //This subfunction accepts the following parameters, in order. | ||

42 | // | ||

43 | // (a) The number to factor. | ||

44 | // | ||

45 | // (b) The number of Miller-Rabin rounds to use when looking | ||

46 | // at the probable primality of the original argument or of | ||

47 | // residues. | ||

48 | // | ||

49 | // (c) The maximum time that should be allowed to elapse before the | ||

50 | // program terminates (perhaps without finding factors). | ||

51 | // This constraint is to protect the server capacity. (In the long | ||

52 | // term, however, the web page may need to be removed--if too many | ||

53 | // people use it, it will bring any server to a crawl.) | ||

54 | // | ||

55 | //OUTPUT RESULTS | ||

56 | //-------------- | ||

57 | //The notation below gives the outputs of the program. In some cases, [i] notation | ||

58 | //is used to indicate line numbers. | ||

59 | // | ||

60 | //[01] An overall success or failure code for the operation. Valid responses | ||

61 | // are: | ||

62 | // S : Success (for this subfunction, the only possible outcome). | ||

63 | // | ||

64 | //[02] The fully normalized first integer entered. This means it has been | ||

65 | // stripped of all weird characters, etc., and also perhaps assigned | ||

66 | // a default value if the original parameter wasn't acceptable. | ||

67 | // This can be used by the PHP script to repopulate form boxes. | ||

68 | //[03] The fully normalized second integer entered. | ||

69 | //[04] The fully normalized second integer entered. | ||

70 | //[05] The type of factor specified on the next line. | ||

71 | // Possibilities are: | ||

72 | // "P" : The factor is definitely prime. | ||

73 | // "C" : The factor is definitely composite. | ||

74 | // "p" : The factor is probably prime, established by | ||

75 | // Miller-Rabin. | ||

76 | // | ||

77 | // Note that the "C" code can only appear if the utility has to give | ||

78 | // up because it runs out of allowed time. This means that the input | ||

79 | // number has not been fully factored or perhaps not factored at all. | ||

80 | // The "C" code can only appear for the last factor or the only factor. | ||

81 | // | ||

82 | // The "p" code can also only appear for the last factor. This occurs when | ||

83 | // either the original input argument is prime as established by | ||

84 | // Miller-Rabin or else a division leaves a result that is similarly | ||

85 | // established by Miller-Rabin. | ||

86 | //[06] The factor itself. | ||

87 | //[07] Its multiplicity. | ||

88 | //[..] Lines 5-7 are repeated for as many factors as are located. | ||

89 | //[..] The second-to-last line will contain "X". | ||

90 | //[..] The last line will contain an "S". | ||

91 | // | ||

92 | //Note the following: | ||

93 | // (a) Any valid output will contain at least 9 lines. | ||

94 | // (b) Any valid output will have a number of lines divisible by 3. | ||

95 | // (c) The number of factors found is (nlines - 6)/3 or alternately | ||

96 | // nlines/3 - 2. | ||

97 | // (d) What happened can be determined using the number of lines plus | ||

98 | // the code of the last factor. | ||

99 | // | ||

100 | //The return value (exit code) from this subfunction is always 0. | ||

101 | // | ||

102 | |||

103 | #define MODULE_SUBFUNC_PFACT_18 | ||

104 | |||

105 | #include <assert.h> | ||

106 | #include <ctype.h> | ||

107 | #include <stddef.h> | ||

108 | #include <stdio.h> | ||

109 | #include <stdlib.h> | ||

110 | #include <string.h> | ||

111 | #include <time.h> | ||

112 | |||

113 | #include <gmp.h> | ||

114 | |||

115 | #include "auxfuncs.h" | ||

116 | #include "subfunc_pfact_18.h" | ||

117 | #include "sieve_eratosthenes.h" | ||

118 | |||

119 | |||

120 | int SUBFUNC_PFACT_18_main(int argc, char *argv[]) | ||

121 | { | ||

122 | //The time snapshot against which we compare to see if we're over | ||

123 | //time budget. | ||

124 | time_t time_snapshot; | ||

125 | |||

126 | //Normalized first and second parameters (the integers). | ||

127 | char *arg1; | ||

128 | char *arg2; | ||

129 | char *arg3; | ||

130 | |||

131 | //Number of Miller-Rabin iterations to establish probable primality. | ||

132 | int miller_rabin_iterations; | ||

133 | |||

134 | //Maximum elapsed time allowed. | ||

135 | int max_time; | ||

136 | |||

137 | //Temporary iteration integers. | ||

138 | int i; | ||

139 | |||

140 | //Mask counter. Only check for timeout periodically, as calling the OS to get time | ||

141 | //is presumed expensive. | ||

142 | int mask_counter; | ||

143 | |||

144 | //The current trial divisor. | ||

145 | unsigned long trial_divisor, new_trial_divisor; | ||

146 | |||

147 | //The current sieve table index. | ||

148 | int sieve_table_index; | ||

149 | |||

150 | //An exit flag kept to remember if we should bail the loop. | ||

151 | int exit_flag; | ||

152 | |||

153 | //The number to factor. | ||

154 | mpz_t number_to_factor; | ||

155 | |||

156 | //The square root limit. We only need to go that far. | ||

157 | mpz_t square_root_limit; | ||

158 | |||

159 | //The multiplicity of any factors we find. | ||

160 | int multiplicity; | ||

161 | |||

162 | //Initialize all of the GMP variables. | ||

163 | mpz_init(number_to_factor); | ||

164 | mpz_init(square_root_limit); | ||

165 | |||

166 | //Grab a time snapshot. | ||

167 | time_snapshot = time(NULL); | ||

168 | |||

169 | //There must be an acceptable number of command-line arguments. If not, | ||

170 | //abort the progam with phony data. | ||

171 | if (argc != 5) | ||

172 | { | ||

173 | printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); | ||

174 | return(0); | ||

175 | } | ||

176 | |||

177 | //Copy the command-line arguments to a safe place where we can manipulate them. | ||

178 | //Leave 2 characters of space in case we assign a "0". | ||

179 | arg1 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[2])) + 1) * sizeof(char)); | ||

180 | arg2 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[3])) + 1) * sizeof(char)); | ||

181 | arg3 = (char *)malloc((AUXFUNCS_size_t_max(1, strlen(argv[4])) + 1) * sizeof(char)); | ||

182 | if ((arg1 == NULL) || (arg2 == NULL) || (arg3 == NULL)) | ||

183 | { | ||

184 | printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); | ||

185 | return(0); | ||

186 | } | ||

187 | strcpy(arg1, argv[2]); | ||

188 | strcpy(arg2, argv[3]); | ||

189 | strcpy(arg3, argv[4]); | ||

190 | |||

191 | //Strip all of the non-digit trash out of the arguments. | ||

192 | AUXFUNCS_remove_non_digits(arg1); | ||

193 | AUXFUNCS_remove_non_digits(arg2); | ||

194 | AUXFUNCS_remove_non_digits(arg3); | ||

195 | |||

196 | //Remove all leading zeros from arguments. | ||

197 | AUXFUNCS_remove_leading_zeros(arg1); | ||

198 | AUXFUNCS_remove_leading_zeros(arg2); | ||

199 | AUXFUNCS_remove_leading_zeros(arg3); | ||

200 | |||

201 | //If an argument is zero length, fill it in with 0. | ||

202 | if (!strlen(arg1)) | ||

203 | strcpy(arg1, "0"); | ||

204 | if (!strlen(arg2)) | ||

205 | strcpy(arg2, "0"); | ||

206 | if (!strlen(arg3)) | ||

207 | strcpy(arg3, "0"); | ||

208 | |||

209 | //We are not allowed to have 0's in this function, so | ||

210 | //abort on zeros. Also, we can't have 1 for a number to check, | ||

211 | //as 1 can't be factored. | ||

212 | if ((!strcmp(arg1, "0")) || (!strcmp(arg2, "0")) || (!strcmp(arg3, "0")) || (!strcmp(arg1, "1"))) | ||

213 | { | ||

214 | printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); | ||

215 | return(0); | ||

216 | } | ||

217 | |||

218 | //If the number to factor exceeds 18 digits, abort. Anything that has | ||

219 | //18 or fewer digits is allowed. | ||

220 | if (strlen(arg1) >18) | ||

221 | { | ||

222 | printf("S\n2\n25\n2\nC\n2\n2\nX\nS\n"); | ||

223 | return(0); | ||

224 | } | ||

225 | |||

226 | //Assign the number to factor. This is definitely a | ||

227 | //GMP type. | ||

228 | mpz_set_str(number_to_factor, arg1, 10); | ||

229 | |||

230 | //Assign the number of Miller-Rabin repetitons to use, subject to an | ||

231 | //absolute maximum of 1000 and minimum of 1. | ||

232 | miller_rabin_iterations = 25; | ||

233 | if (strlen(arg2) > 3) | ||

234 | { | ||

235 | miller_rabin_iterations = 1000; | ||

236 | } | ||

237 | else | ||

238 | { | ||

239 | sscanf(arg2, "%d", &miller_rabin_iterations); | ||

240 | if (miller_rabin_iterations < 1) | ||

241 | miller_rabin_iterations = 1; | ||

242 | } | ||

243 | |||

244 | //Assign the maximum time allowed, subject to a maximum of 1000 and minimum of 3. | ||

245 | if (strlen(arg3) > 3) | ||

246 | { | ||

247 | max_time = 1000; | ||

248 | } | ||

249 | else | ||

250 | { | ||

251 | sscanf(arg3, "%d", &max_time); | ||

252 | if (max_time < 3) | ||

253 | max_time = 3; | ||

254 | } | ||

255 | |||

256 | //Output the header information before beginning the search. | ||

257 | printf("S\n"); | ||

258 | mpz_out_str(stdout, 10, number_to_factor); | ||

259 | printf("\n"); | ||

260 | printf("%d\n", miller_rabin_iterations); | ||

261 | printf("%d\n", max_time); | ||

262 | |||

263 | //Loop through the list of divisors we should try before starting the sieve. | ||

264 | //This won't possibly exceed our time, so do it. | ||

265 | for (i=0; i<SIEVE_ERATOSTHENES_N_SIEVE_FACTORS; i++) | ||

266 | { | ||

267 | multiplicity = 0; | ||

268 | |||

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

270 | |||

271 | //Factor out all occurrences that we can. | ||

272 | while (mpz_divisible_ui_p(number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i])) | ||

273 | { | ||

274 | mpz_divexact_ui(number_to_factor, number_to_factor, SIEVE_ERATOSTHENES_sieve_factors[i]); | ||

275 | multiplicity++; | ||

276 | } | ||

277 | |||

278 | //Now output the information record to the stdout if we could factor it. | ||

279 | if (multiplicity) | ||

280 | { | ||

281 | printf("P\n"); | ||

282 | printf("%u\n", SIEVE_ERATOSTHENES_sieve_factors[i]); | ||

283 | printf("%d\n", multiplicity); | ||

284 | } | ||

285 | } | ||

286 | |||

287 | //Gear up for tabulated operation as we sieve. | ||

288 | new_trial_divisor = trial_divisor = SIEVE_ERATOSTHENES_FIRST_TRIAL_DIVISOR; | ||

289 | sieve_table_index = SIEVE_ERATOSTHENES_FIRST_SIEVE_INDEX; | ||

290 | mask_counter = 0; | ||

291 | |||

292 | //Set the exit flag initially. One thing that could have | ||

293 | //caused us to be ready to exit already is if we brought the | ||

294 | //number to factor down to 1. We check this now because | ||

295 | //otherwise only check it when have factored something out and | ||

296 | //seive loop would run indefinitely if didn't do this. | ||

297 | exit_flag = 0; | ||

298 | if (!mpz_cmp_ui(number_to_factor, 1)) | ||

299 | exit_flag = 1; | ||

300 | |||

301 | //Effectively, we may have broken down the number to factor | ||

302 | //with just a few operations. So, effectively, we are | ||

303 | //starting afresh here. We only need to proceed to the | ||

304 | //square root of the current number to factor (not the | ||

305 | //original one. | ||

306 | mpz_sqrt(square_root_limit, number_to_factor); | ||

307 | |||

308 | //If the last value from the numbers to try before we start | ||

309 | //the sieve is already past the square root limit, then whatever | ||

310 | //remains is definitely prime, and we can exit. | ||

311 | if (mpz_cmp_ui(square_root_limit, SIEVE_ERATOSTHENES_sieve_factors[SIEVE_ERATOSTHENES_N_SIEVE_FACTORS - 1]) <= 0) | ||

312 | { | ||

313 | if (mpz_cmp_ui(number_to_factor, 1)) | ||

314 | { | ||

315 | exit_flag = 1; | ||

316 | printf("P\n"); | ||

317 | mpz_out_str(stdout, 10, number_to_factor); | ||

318 | printf("\n"); | ||

319 | printf("1\n"); | ||

320 | goto done; | ||

321 | } | ||

322 | } | ||

323 | |||

324 | //If Miller-Rabin says that the remaining number is probably prime, that is good | ||

325 | //enough. | ||

326 | if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1) | ||

327 | { | ||

328 | exit_flag = 1; | ||

329 | printf("p\n"); | ||

330 | mpz_out_str(stdout, 10, number_to_factor); | ||

331 | printf("\n"); | ||

332 | printf("1\n"); | ||

333 | goto done; | ||

334 | } | ||

335 | |||

336 | //Loop entry condition. | ||

337 | multiplicity = 0; | ||

338 | |||

339 | //This is the main loop. We only check for elapsed time | ||

340 | //infrequently. | ||

341 | while (!exit_flag) | ||

342 | { | ||

343 | //Replace the trial divisor. | ||

344 | trial_divisor = new_trial_divisor; | ||

345 | |||

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

347 | |||

348 | //Factor out all occurrences that we can of the trial divisor, if we can. | ||

349 | while (mpz_divisible_ui_p(number_to_factor, trial_divisor)) | ||

350 | { | ||

351 | mpz_divexact_ui(number_to_factor, number_to_factor, trial_divisor); | ||

352 | multiplicity++; | ||

353 | } | ||

354 | |||

355 | //Now output the information record to the stdout if we could factor it. | ||

356 | if (multiplicity) | ||

357 | { | ||

358 | printf("P\n"); | ||

359 | printf("%u\n", trial_divisor); | ||

360 | printf("%d\n", multiplicity); | ||

361 | |||

362 | //We need to calculate a new square root bound. | ||

363 | mpz_sqrt(square_root_limit, number_to_factor); | ||

364 | |||

365 | //If we are over the square root bound, the remaining number to factor | ||

366 | //is prime and we should exit. | ||

367 | if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) | ||

368 | exit_flag = 1; | ||

369 | |||

370 | //Are we over the time bound? We might has well check here, | ||

371 | //because finding prime factors is rare. | ||

372 | if ((time(NULL) - time_snapshot) > max_time) | ||

373 | exit_flag = 1; | ||

374 | |||

375 | //Does the Miller-Rabin test say that the remaining number | ||

376 | //is with near perfect certainty prime? | ||

377 | if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) >= 1) | ||

378 | exit_flag = 1; | ||

379 | |||

380 | //Is the remaining number "1", indicating we've factored it fully? | ||

381 | if (!mpz_cmp_ui(number_to_factor, 1)) | ||

382 | exit_flag = 1; | ||

383 | |||

384 | //Set multiplicity to zero again so don't re-enter until successful | ||

385 | //division again. | ||

386 | multiplicity = 0; | ||

387 | } | ||

388 | |||

389 | //Advance to our next trial divisor. | ||

390 | new_trial_divisor = trial_divisor + SIEVE_ERATOSTHENES_sieve[sieve_table_index]; | ||

391 | |||

392 | //If the new is < the old, we've rolled over. This means an exit is necessary. | ||

393 | if (new_trial_divisor < trial_divisor) | ||

394 | exit_flag = 1; | ||

395 | |||

396 | //Advance the sieve index. | ||

397 | sieve_table_index = (sieve_table_index + 1) % SIEVE_ERATOSTHENES_N_SIEVE; | ||

398 | |||

399 | //This is our only chance to check for termination conditions that | ||

400 | //don't come about from a successful division. But we don't do | ||

401 | //this often. | ||

402 | mask_counter++; | ||

403 | if (!(mask_counter & 0xFFFFF)) | ||

404 | { | ||

405 | //First, have we exceeded the square root bound? | ||

406 | if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) | ||

407 | exit_flag = 1; | ||

408 | //Second, are we over the time budget? | ||

409 | if ((time(NULL) - time_snapshot) > max_time) | ||

410 | exit_flag = 1; | ||

411 | } | ||

412 | } | ||

413 | |||

414 | //If we've made it out of the loop, there are a variety of reasons for that. | ||

415 | //Find the right one and close up. | ||

416 | if (!mpz_cmp_ui(number_to_factor, 1)) | ||

417 | { | ||

418 | //We divided the number down to 1. There is nothing further to do. | ||

419 | } | ||

420 | else if (mpz_cmp_ui(square_root_limit, trial_divisor) <= 0) | ||

421 | { | ||

422 | //We are at or over the square root limit. The remaining number is definitely prime. | ||

423 | printf("P\n"); | ||

424 | mpz_out_str(stdout, 10, number_to_factor); | ||

425 | printf("\n"); | ||

426 | printf("1\n"); | ||

427 | } | ||

428 | else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 0) | ||

429 | { | ||

430 | //Miller-Rabin says the remaining number is definitely composite. | ||

431 | printf("C\n"); | ||

432 | mpz_out_str(stdout, 10, number_to_factor); | ||

433 | printf("\n"); | ||

434 | printf("1\n"); | ||

435 | } | ||

436 | else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 1) | ||

437 | { | ||

438 | //Miller-Rabin says the remaining number is probably prime. | ||

439 | printf("p\n"); | ||

440 | mpz_out_str(stdout, 10, number_to_factor); | ||

441 | printf("\n"); | ||

442 | printf("1\n"); | ||

443 | } | ||

444 | else if (mpz_probab_prime_p(number_to_factor, miller_rabin_iterations) == 2) | ||

445 | { | ||

446 | //Miller-Rabin says the remaining number is definitely prime. | ||

447 | printf("P\n"); | ||

448 | mpz_out_str(stdout, 10, number_to_factor); | ||

449 | printf("\n"); | ||

450 | printf("1\n"); | ||

451 | } | ||

452 | |||

453 | //Output the invariant footer information. | ||

454 | done: | ||

455 | printf("X\n"); | ||

456 | printf("S\n"); | ||

457 | |||

458 | //Always return 0. | ||

459 | return(0); | ||

460 | } | ||

461 | |||

462 | //******************************************************************************** | ||

463 | // End of SUBFUNC_PFACT_18.C. | ||

464 | //******************************************************************************** |

Name | Value |
---|---|

svn:eol-style |
native |

svn:keywords |
Header |

dashley@gmail.com | ViewVC Help |

Powered by ViewVC 1.1.25 |