source/libs/gmp/gmp-src/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-2002, 2005, 2014 Free Software
  13 Foundation, Inc.
  14
  15 Contributed by John Amanatides.
  16
  17 This file is part of the GNU MP Library.
  18
  19 The GNU MP Library is free software; you can redistribute it and/or modify
  20 it under the terms of either:
  21
  22   * the GNU Lesser General Public License as published by the Free
  23     Software Foundation; either version 3 of the License, or (at your
  24     option) any later version.
  25
  26 or
  27
  28   * the GNU General Public License as published by the Free Software
  29     Foundation; either version 2 of the License, or (at your option) any
  30     later version.
  31
  32 or both in parallel, as here.
  33
  34 The GNU MP Library is distributed in the hope that it will be useful, but
  35 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  36 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  37 for more details.
  38
  39 You should have received copies of the GNU General Public License and the
  40 GNU Lesser General Public License along with the GNU MP Library.  If not,
  41 see https://www.gnu.org/licenses/.  */
  42
  43 #include "gmp.h"
  44 #include "gmp-impl.h"
  45
  46 static int millerrabin (mpz_srcptr, mpz_srcptr,
  47                         mpz_ptr, mpz_ptr,
  48                         mpz_srcptr, unsigned long int);
  49
  50 int
  51 mpz_millerrabin (mpz_srcptr n, int reps)
  52 {
  53   int r;
  54   mpz_t nm1, nm3, x, y, q;
  55   unsigned long int k;
  56   gmp_randstate_t rstate;
  57   int is_prime;
  58   TMP_DECL;
  59   TMP_MARK;
  60
  61   MPZ_TMP_INIT (nm1, SIZ (n) + 1);
  62   mpz_sub_ui (nm1, n, 1L);
  63
  64   MPZ_TMP_INIT (x, SIZ (n) + 1);
  65   MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
  66
  67   /* Perform a Fermat test.  */
  68   mpz_set_ui (x, 210L);
  69   mpz_powm (y, x, nm1, n);
  70   if (mpz_cmp_ui (y, 1L) != 0)
  71     {
  72       TMP_FREE;
  73       return 0;
  74     }
  75
  76   MPZ_TMP_INIT (q, SIZ (n));
  77
  78   /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
  79   k = mpz_scan1 (nm1, 0L);
  80   mpz_tdiv_q_2exp (q, nm1, k);
  81
  82   /* n-3 */
  83   MPZ_TMP_INIT (nm3, SIZ (n) + 1);
  84   mpz_sub_ui (nm3, n, 3L);
  85   ASSERT (mpz_cmp_ui (nm3, 1L) >= 0);
  86
  87   gmp_randinit_default (rstate);
  88
  89   is_prime = 1;
  90   for (r = 0; r < reps && is_prime; r++)
  91     {
  92       /* 2 to n-2 inclusive, don't want 1, 0 or -1 */
  93       mpz_urandomm (x, rstate, nm3);
  94       mpz_add_ui (x, x, 2L);
  95
  96       is_prime = millerrabin (n, nm1, x, y, q, k);
  97     }
  98
  99   gmp_randclear (rstate);
 100
 101   TMP_FREE;
 102   return is_prime;
 103 }
 104
 105 static int
 106 millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
 107              mpz_srcptr q, unsigned long int k)
 108 {
 109   unsigned long int i;
 110
 111   mpz_powm (y, x, q, n);
 112
 113   if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
 114     return 1;
 115
 116   for (i = 1; i < k; i++)
 117     {
 118       mpz_powm_ui (y, y, 2L, n);
 119       if (mpz_cmp (y, nm1) == 0)
 120         return 1;
 121       /* y == 1 means that the previous y was a non-trivial square root
 122          of 1 (mod n). y == 0 means that n is a power of the base.
 123          In either case, n is not prime. */
 124       if (mpz_cmp_ui (y, 1L) <= 0)
 125         return 0;
 126     }
 127   return 0;
 128 }