games/primes/spsp.c

   1 /*      $NetBSD: spsp.c,v 1.1 2014/10/02 21:36:37 ast 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.1 2014/10/02 21:36:37 ast 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 <= a, b < 2^63, 0 < n < 2^63. */
  50 static uint64_t
  51 mulmod(uint64_t a, uint64_t b, uint64_t n)
  52 {
  53         uint64_t x = 0;
  54
  55         while (b != 0) {
  56                 if (b & 1)
  57                         x = (x + a) % n;
  58                 a = (a + a) % n;
  59                 b >>= 1;
  60         }
  61
  62         return (x);
  63 }
  64
  65 /* Return a^r % n, where 0 <= a < 2^63, 0 < n < 2^63. */
  66 static uint64_t
  67 powmod(uint64_t a, uint64_t r, uint64_t n)
  68 {
  69         uint64_t x = 1;
  70
  71         while (r != 0) {
  72                 if (r & 1)
  73                         x = mulmod(a, x, n);
  74                 a = mulmod(a, a, n);
  75                 r >>= 1;
  76         }
  77
  78         return (x);
  79 }
  80
  81 /* Return non-zero if n is a strong pseudoprime to base p. */
  82 static int
  83 spsp(uint64_t n, uint64_t p)
  84 {
  85         uint64_t x;
  86         uint64_t r = n - 1;
  87         int k = 0;
  88
  89         /* Compute n - 1 = 2^k * r. */
  90         while ((r & 1) == 0) {
  91                 k++;
  92                 r >>= 1;
  93         }
  94
  95         /* Compute x = p^r mod n.  If x = 1, n is a p-spsp. */
  96         x = powmod(p, r, n);
  97         if (x == 1)
  98                 return (1);
  99
 100         /* Compute x^(2^i) for 0 <= i < n.  If any are -1, n is a p-spsp. */
 101         while (k > 0) {
 102                 if (x == n - 1)
 103                         return (1);
 104                 x = powmod(x, 2, n);
 105                 k--;
 106         }
 107
 108         /* Not a p-spsp. */
 109         return (0);
 110 }
 111
 112 /* Test for primality using strong pseudoprime tests. */
 113 int
 114 isprime(uint64_t _n)
 115 {
 116         uint64_t n = _n;
 117
 118         /*
 119          * Values from:
 120          * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
 121          * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
 122          */
 123
 124         /* No SPSPs to base 2 less than 2047. */
 125         if (!spsp(n, 2))
 126                 return (0);
 127         if (n < 2047ULL)
 128                 return (1);
 129
 130         /* No SPSPs to bases 2,3 less than 1373653. */
 131         if (!spsp(n, 3))
 132                 return (0);
 133         if (n < 1373653ULL)
 134                 return (1);
 135
 136         /* No SPSPs to bases 2,3,5 less than 25326001. */
 137         if (!spsp(n, 5))
 138                 return (0);
 139         if (n < 25326001ULL)
 140                 return (1);
 141
 142         /* No SPSPs to bases 2,3,5,7 less than 3215031751. */
 143         if (!spsp(n, 7))
 144                 return (0);
 145         if (n < 3215031751ULL)
 146                 return (1);
 147
 148         /*
 149          * Values from:
 150          * G. Jaeschke, On strong pseudoprimes to several bases,
 151          * Math. Comp. 61(204):915-926, 1993.
 152          */
 153
 154         /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
 155         if (!spsp(n, 11))
 156                 return (0);
 157         if (n < 2152302898747ULL)
 158                 return (1);
 159
 160         /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
 161         if (!spsp(n, 13))
 162                 return (0);
 163         if (n < 3474749660383ULL)
 164                 return (1);
 165
 166         /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
 167         if (!spsp(n, 17))
 168                 return (0);
 169         if (n < 341550071728321ULL)
 170                 return (1);
 171
 172         /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
 173         if (!spsp(n, 19))
 174                 return (0);
 175         if (n < 341550071728321ULL)
 176                 return (1);
 177
 178         /*
 179          * Value from:
 180          * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
 181          * bases, Math. Comp. 83(290):2915-2924, 2014.
 182          */
 183
 184         /* No SPSPs to bases 2..23 less than 3825123056546413051. */
 185         if (!spsp(n, 23))
 186                 return (0);
 187         if (n < 3825123056546413051)
 188                 return (1);
 189
 190         /* We can't handle values larger than this. */
 191         assert(n <= SPSPMAX);
 192
 193         /* UNREACHABLE */
 194         return (0);
 195 }