usr.bin/primes/spsp.c

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