#
source:
trunk/third/gmp/mpz/pprime_p.c
@
18191

Revision 18191, 4.0 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* mpz_probab_prime_p -- |

2 | An implementation of the probabilistic primality test found in Knuth's |

3 | Seminumerical Algorithms book. If the function mpz_probab_prime_p() |

4 | returns 0 then n is not prime. If it returns 1, then n is 'probably' |

5 | prime. If it returns 2, n is surely prime. The probability of a false |

6 | positive is (1/4)**reps, where reps is the number of internal passes of the |

7 | probabilistic algorithm. Knuth indicates that 25 passes are reasonable. |

8 | |

9 | Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software |

10 | Foundation, Inc. Miller-Rabin code contributed by John Amanatides. |

11 | |

12 | This file is part of the GNU MP Library. |

13 | |

14 | The GNU MP Library is free software; you can redistribute it and/or modify |

15 | it under the terms of the GNU Lesser General Public License as published by |

16 | the Free Software Foundation; either version 2.1 of the License, or (at your |

17 | option) any later version. |

18 | |

19 | The GNU MP Library is distributed in the hope that it will be useful, but |

20 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

21 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

22 | License for more details. |

23 | |

24 | You should have received a copy of the GNU Lesser General Public License |

25 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

26 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

27 | MA 02111-1307, USA. */ |

28 | |

29 | #include "gmp.h" |

30 | #include "gmp-impl.h" |

31 | #include "longlong.h" |

32 | |

33 | static int isprime _PROTO ((unsigned long int t)); |

34 | |

35 | int |

36 | mpz_probab_prime_p (mpz_srcptr n, int reps) |

37 | { |

38 | mp_limb_t r; |

39 | |

40 | /* Handle small and negative n. */ |

41 | if (mpz_cmp_ui (n, 1000000L) <= 0) |

42 | { |

43 | int is_prime; |

44 | if (mpz_sgn (n) < 0) |

45 | { |

46 | /* Negative number. Negate and call ourselves. */ |

47 | mpz_t n2; |

48 | mpz_init (n2); |

49 | mpz_neg (n2, n); |

50 | is_prime = mpz_probab_prime_p (n2, reps); |

51 | mpz_clear (n2); |

52 | return is_prime; |

53 | } |

54 | is_prime = isprime (mpz_get_ui (n)); |

55 | return is_prime ? 2 : 0; |

56 | } |

57 | |

58 | /* If n is now even, it is not a prime. */ |

59 | if ((mpz_get_ui (n) & 1) == 0) |

60 | return 0; |

61 | |

62 | #if defined (PP) |

63 | /* Check if n has small factors. */ |

64 | #if defined (PP_INVERTED) |

65 | r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP, |

66 | (mp_limb_t) PP_INVERTED); |

67 | #else |

68 | r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); |

69 | #endif |

70 | if (r % 3 == 0 |

71 | #if BITS_PER_MP_LIMB >= 4 |

72 | || r % 5 == 0 |

73 | #endif |

74 | #if BITS_PER_MP_LIMB >= 8 |

75 | || r % 7 == 0 |

76 | #endif |

77 | #if BITS_PER_MP_LIMB >= 16 |

78 | || r % 11 == 0 || r % 13 == 0 |

79 | #endif |

80 | #if BITS_PER_MP_LIMB >= 32 |

81 | || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 |

82 | #endif |

83 | #if BITS_PER_MP_LIMB >= 64 |

84 | || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0 |

85 | || r % 47 == 0 || r % 53 == 0 |

86 | #endif |

87 | ) |

88 | { |

89 | return 0; |

90 | } |

91 | #endif /* PP */ |

92 | |

93 | /* Do more dividing. We collect small primes, using umul_ppmm, until we |

94 | overflow a single limb. We divide our number by the small primes product, |

95 | and look for factors in the remainder. */ |

96 | { |

97 | unsigned long int ln2; |

98 | unsigned long int q; |

99 | mp_limb_t p1, p0, p; |

100 | unsigned int primes[15]; |

101 | int nprimes; |

102 | |

103 | nprimes = 0; |

104 | p = 1; |

105 | ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; |

106 | for (q = PP_FIRST_OMITTED; q < ln2; q += 2) |

107 | { |

108 | if (isprime (q)) |

109 | { |

110 | umul_ppmm (p1, p0, p, q); |

111 | if (p1 != 0) |

112 | { |

113 | r = mpn_mod_1 (PTR(n), SIZ(n), p); |

114 | while (--nprimes >= 0) |

115 | if (r % primes[nprimes] == 0) |

116 | { |

117 | ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0); |

118 | return 0; |

119 | } |

120 | p = q; |

121 | nprimes = 0; |

122 | } |

123 | else |

124 | { |

125 | p = p0; |

126 | } |

127 | primes[nprimes++] = q; |

128 | } |

129 | } |

130 | } |

131 | |

132 | /* Perform a number of Miller-Rabin tests. */ |

133 | return mpz_millerrabin (n, reps); |

134 | } |

135 | |

136 | static int |

137 | isprime (unsigned long int t) |

138 | { |

139 | unsigned long int q, r, d; |

140 | |

141 | if (t < 3 || (t & 1) == 0) |

142 | return t == 2; |

143 | |

144 | for (d = 3, r = 1; r != 0; d += 2) |

145 | { |

146 | q = t / d; |

147 | r = t - q * d; |

148 | if (q < d) |

149 | return 1; |

150 | } |

151 | return 0; |

152 | } |

**Note:**See TracBrowser for help on using the repository browser.