contrib/gmp/mpz/millerrabin.c

   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, 2002, 2005 Free
  13 Software 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 3 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.  If not, see http://www.gnu.org/licenses/.  */
  29
  30 #include "gmp.h"
  31 #include "gmp-impl.h"
  32
  33 static int millerrabin __GMP_PROTO ((mpz_srcptr, mpz_srcptr,
  34                                      mpz_ptr, mpz_ptr,
  35                                      mpz_srcptr, unsigned long int));
  36
  37 int
  38 mpz_millerrabin (mpz_srcptr n, int reps)
  39 {
  40   int r;
  41   mpz_t nm1, nm3, x, y, q;
  42   unsigned long int k;
  43   gmp_randstate_t rstate;
  44   int is_prime;
  45   TMP_DECL;
  46   TMP_MARK;
  47
  48   MPZ_TMP_INIT (nm1, SIZ (n) + 1);
  49   mpz_sub_ui (nm1, n, 1L);
  50
  51   MPZ_TMP_INIT (x, SIZ (n) + 1);
  52   MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
  53
  54   /* Perform a Fermat test.  */
  55   mpz_set_ui (x, 210L);
  56   mpz_powm (y, x, nm1, n);
  57   if (mpz_cmp_ui (y, 1L) != 0)
  58     {
  59       TMP_FREE;
  60       return 0;
  61     }
  62
  63   MPZ_TMP_INIT (q, SIZ (n));
  64
  65   /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
  66   k = mpz_scan1 (nm1, 0L);
  67   mpz_tdiv_q_2exp (q, nm1, k);
  68
  69   /* n-3 */
  70   MPZ_TMP_INIT (nm3, SIZ (n) + 1);
  71   mpz_sub_ui (nm3, n, 3L);
  72   ASSERT (mpz_cmp_ui (nm3, 1L) >= 0);
  73
  74   gmp_randinit_default (rstate);
  75
  76   is_prime = 1;
  77   for (r = 0; r < reps && is_prime; r++)
  78     {
  79       /* 2 to n-2 inclusive, don't want 1, 0 or -1 */
  80       mpz_urandomm (x, rstate, nm3);
  81       mpz_add_ui (x, x, 2L);
  82
  83       is_prime = millerrabin (n, nm1, x, y, q, k);
  84     }
  85
  86   gmp_randclear (rstate);
  87
  88   TMP_FREE;
  89   return is_prime;
  90 }
  91
  92 static int
  93 millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
  94              mpz_srcptr q, unsigned long int k)
  95 {
  96   unsigned long int i;
  97
  98   mpz_powm (y, x, q, n);
  99
 100   if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
 101     return 1;
 102
 103   for (i = 1; i < k; i++)
 104     {
 105       mpz_powm_ui (y, y, 2L, n);
 106       if (mpz_cmp (y, nm1) == 0)
 107         return 1;
 108       if (mpz_cmp_ui (y, 1L) == 0)
 109         return 0;
 110     }
 111   return 0;
 112 }