primes/spsp.c

   1 /*      $NetBSD: spsp.c,v 1.2 2018/02/03 15:40:29 christos Exp $        */
   2
   3 /*-
   4  * Copyright (c) 2014 Colin Percival
   5  * All rights reserved.
   6  *
   7  * Redistribution and use in source and binary forms, with or without
   8  * modification, are permitted provided that the following conditions
   9  * are met:
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  * 2. Redistributions in binary form must reproduce the above copyright
  13  *    notice, this list of conditions and the following disclaimer in the
  14  *    documentation and/or other materials provided with the distribution.
  15  *
  16  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  20  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  21  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  22  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  23  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  24  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  25  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  26  * SUCH DAMAGE.
  27  */
  28
  29 #include <sys/cdefs.h>
  30 #ifndef lint
  31 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
  32  The Regents of the University of California.  All rights reserved.");
  33 #endif /* not lint */
  34
  35 #ifndef lint
  36 #if 0
  37 static char sccsid[] = "@(#)primes.c    8.5 (Berkeley) 5/10/95";
  38 #else
  39 __RCSID("$NetBSD: spsp.c,v 1.2 2018/02/03 15:40:29 christos Exp $");
  40 #endif
  41 #endif /* not lint */
  42
  43 #include <assert.h>
  44 #include <stddef.h>
  45 #include <stdint.h>
  46
  47 #include "primes.h"
  48
  49 /* Return a * b % n, where 0 <= n. */
  50 static uint64_t
  51 mulmod(uint64_t a, uint64_t b, uint64_t n)
  52 {
  53         uint64_t x = 0;
  54         uint64_t an = a % n;
  55
  56         while (b != 0) {
  57                 if (b & 1) {
  58                         x += an;
  59                         if ((x < an) || (x >= n))
  60                                 x -= n;
  61                 }
  62                 if (an + an < an)
  63                         an = an + an - n;
  64                 else if (an + an >= n)
  65                         an = an + an - n;
  66                 else
  67                         an = an + an;
  68
  69                 b >>= 1;
  70         }
  71
  72         return (x);
  73 }
  74
  75 /* Return a^r % n, where 0 < n. */
  76 static uint64_t
  77 powmod(uint64_t a, uint64_t r, uint64_t n)
  78 {
  79         uint64_t x = 1;
  80
  81         while (r != 0) {
  82                 if (r & 1)
  83                         x = mulmod(a, x, n);
  84                 a = mulmod(a, a, n);
  85                 r >>= 1;
  86         }
  87
  88         return (x);
  89 }
  90
  91 /* Return non-zero if n is a strong pseudoprime to base p. */
  92 static int
  93 spsp(uint64_t n, uint64_t p)
  94 {
  95         uint64_t x;
  96         uint64_t r = n - 1;
  97         int k = 0;
  98
  99         /* Compute n - 1 = 2^k * r. */
 100         while ((r & 1) == 0) {
 101                 k++;
 102                 r >>= 1;
 103         }
 104
 105         /* Compute x = p^r mod n.  If x = 1, n is a p-spsp. */
 106         x = powmod(p, r, n);
 107         if (x == 1)
 108                 return (1);
 109
 110         /* Compute x^(2^i) for 0 <= i < n.  If any are -1, n is a p-spsp. */
 111         while (k > 0) {
 112                 if (x == n - 1)
 113                         return (1);
 114                 x = powmod(x, 2, n);
 115                 k--;
 116         }
 117
 118         /* Not a p-spsp. */
 119         return (0);
 120 }
 121
 122 /* Test for primality using strong pseudoprime tests. */
 123 int
 124 isprime(uint64_t _n)
 125 {
 126         uint64_t n = _n;
 127
 128         /*
 129          * Values from:
 130          * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
 131          * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
 132          */
 133
 134         /* No SPSPs to base 2 less than 2047. */
 135         if (!spsp(n, 2))
 136                 return (0);
 137         if (n < 2047ULL)
 138                 return (1);
 139
 140         /* No SPSPs to bases 2,3 less than 1373653. */
 141         if (!spsp(n, 3))
 142                 return (0);
 143         if (n < 1373653ULL)
 144                 return (1);
 145
 146         /* No SPSPs to bases 2,3,5 less than 25326001. */
 147         if (!spsp(n, 5))
 148                 return (0);
 149         if (n < 25326001ULL)
 150                 return (1);
 151
 152         /* No SPSPs to bases 2,3,5,7 less than 3215031751. */
 153         if (!spsp(n, 7))
 154                 return (0);
 155         if (n < 3215031751ULL)
 156                 return (1);
 157
 158         /*
 159          * Values from:
 160          * G. Jaeschke, On strong pseudoprimes to several bases,
 161          * Math. Comp. 61(204):915-926, 1993.
 162          */
 163
 164         /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
 165         if (!spsp(n, 11))
 166                 return (0);
 167         if (n < 2152302898747ULL)
 168                 return (1);
 169
 170         /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
 171         if (!spsp(n, 13))
 172                 return (0);
 173         if (n < 3474749660383ULL)
 174                 return (1);
 175
 176         /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
 177         if (!spsp(n, 17))
 178                 return (0);
 179         if (n < 341550071728321ULL)
 180                 return (1);
 181
 182         /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
 183         if (!spsp(n, 19))
 184                 return (0);
 185         if (n < 341550071728321ULL)
 186                 return (1);
 187
 188         /*
 189          * Value from:
 190          * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
 191          * bases, Math. Comp. 83(290):2915-2924, 2014.
 192          */
 193
 194         /* No SPSPs to bases 2..23 less than 3825123056546413051. */
 195         if (!spsp(n, 23))
 196                 return (0);
 197         if (n < 3825123056546413051)
 198                 return (1);
 199         /*
 200          * Value from:
 201          * J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
 202          * bases, Math. Comp. 86(304):985-1003, 2017.
 203          */
 204
 205        /* No SPSPs to bases 2..37 less than 318665857834031151167461. */
 206        if (!spsp(n, 29))
 207                return (0);
 208        if (!spsp(n, 31))
 209                return (0);
 210        if (!spsp(n, 37))
 211                return (0);
 212
 213        /* All 64-bit values are less than 318665857834031151167461. */
 214        return (1);
 215 }