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

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

Line | |
---|---|

1 | /* mpz_millerrabin(n,reps) -- An implementation of the probabilistic primality |

2 | test found in Knuth's Seminumerical Algorithms book. If the function |

3 | mpz_millerrabin() returns 0 then n is not prime. If it returns 1, then n is |

4 | 'probably' prime. The probability of a false positive is (1/4)**reps, where |

5 | reps is the number of internal passes of the probabilistic algorithm. Knuth |

6 | indicates that 25 passes are reasonable. |

7 | |

8 | THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST |

9 | CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN |

10 | FUTURE GNU MP RELEASES. |

11 | |

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

13 | Foundation, Inc. Contributed by John Amanatides. |

14 | |

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

16 | |

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

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

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

20 | option) any later version. |

21 | |

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

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

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

25 | License for more details. |

26 | |

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

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

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

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

31 | |

32 | #include "gmp.h" |

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

34 | |

35 | static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1, |

36 | mpz_ptr x, mpz_ptr y, |

37 | mpz_srcptr q, unsigned long int k)); |

38 | |

39 | int |

40 | mpz_millerrabin (mpz_srcptr n, int reps) |

41 | { |

42 | int r; |

43 | mpz_t nm1, x, y, q; |

44 | unsigned long int k; |

45 | gmp_randstate_t rstate; |

46 | int is_prime; |

47 | TMP_DECL (marker); |

48 | TMP_MARK (marker); |

49 | |

50 | MPZ_TMP_INIT (nm1, SIZ (n) + 1); |

51 | mpz_sub_ui (nm1, n, 1L); |

52 | |

53 | MPZ_TMP_INIT (x, SIZ (n)); |

54 | MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ |

55 | |

56 | /* Perform a Fermat test. */ |

57 | mpz_set_ui (x, 210L); |

58 | mpz_powm (y, x, nm1, n); |

59 | if (mpz_cmp_ui (y, 1L) != 0) |

60 | { |

61 | TMP_FREE (marker); |

62 | return 0; |

63 | } |

64 | |

65 | MPZ_TMP_INIT (q, SIZ (n)); |

66 | |

67 | /* Find q and k, where q is odd and n = 1 + 2**k * q. */ |

68 | k = mpz_scan1 (nm1, 0L); |

69 | mpz_tdiv_q_2exp (q, nm1, k); |

70 | |

71 | gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L); |

72 | |

73 | is_prime = 1; |

74 | for (r = 0; r < reps && is_prime; r++) |

75 | { |

76 | do |

77 | mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1); |

78 | while (mpz_cmp_ui (x, 1L) <= 0); |

79 | |

80 | is_prime = millerrabin (n, nm1, x, y, q, k); |

81 | } |

82 | |

83 | gmp_randclear (rstate); |

84 | |

85 | TMP_FREE (marker); |

86 | return is_prime; |

87 | } |

88 | |

89 | static int |

90 | millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y, |

91 | mpz_srcptr q, unsigned long int k) |

92 | { |

93 | unsigned long int i; |

94 | |

95 | mpz_powm (y, x, q, n); |

96 | |

97 | if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0) |

98 | return 1; |

99 | |

100 | for (i = 1; i < k; i++) |

101 | { |

102 | mpz_powm_ui (y, y, 2L, n); |

103 | if (mpz_cmp (y, nm1) == 0) |

104 | return 1; |

105 | if (mpz_cmp_ui (y, 1L) == 0) |

106 | return 0; |

107 | } |

108 | return 0; |

109 | } |

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