Line data Source code
1 : /* sp.c - "small prime" functions that don't need to be inlined
2 :
3 : Copyright 2005, 2006, 2007, 2008, 2009, 2010 Dave Newman, Jason Papadopoulos,
4 : Alexander Kruppa, Paul Zimmermann.
5 :
6 : The SP Library is free software; you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation; either version 3 of the License, or (at your
9 : option) any later version.
10 :
11 : The SP Library is distributed in the hope that it will be useful, but
12 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 : License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with the SP Library; see the file COPYING.LIB. If not, write to
18 : the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
19 : MA 02110-1301, USA. */
20 :
21 : #include <stdio.h> /* for stderr */
22 : #include <stdlib.h>
23 : #include "sp.h"
24 :
25 : /* Test if m is a base "a" strong probable prime */
26 :
27 : static int
28 689684 : sp_spp (sp_t a, sp_t m, sp_t d)
29 : {
30 : sp_t r, s, t, e;
31 :
32 : /* since this routine is called only with m >= SP_MIN >= 2^30,
33 : and a <= 61, we cannot have m = a */
34 : ASSERT(a < m);
35 :
36 : /* Set e * 2^s = m-1, e odd */
37 6731518 : for (s = 0, e = m - 1; !(e & 1); s++, e >>= 1);
38 :
39 689684 : t = sp_pow (a, e, m, d);
40 :
41 689684 : if (t == 1)
42 7761 : return 1;
43 :
44 6182659 : for (r = 0; r < s; r++)
45 : {
46 5755335 : if (t == m - 1)
47 254599 : return 1;
48 :
49 5500736 : t = sp_sqr (t, m, d);
50 : }
51 :
52 427324 : return 0;
53 : }
54 :
55 : /* Test if x is a prime, return 1 if it is. Note this only works on sp's,
56 : i.e. we need the top bit of x set.
57 : Assume x is odd and x >= SP_MIN.
58 : */
59 : int
60 453560 : sp_prime (sp_t x)
61 : {
62 : sp_t d;
63 :
64 : ASSERT (x & 1);
65 : ASSERT (x >= SP_MIN);
66 :
67 453560 : sp_reciprocal (d, x);
68 :
69 : if (SP_NUMB_BITS <= 32)
70 : {
71 : /* 32-bit primality test
72 : * See http://primes.utm.edu/prove/prove2_3.html */
73 :
74 : if (!sp_spp (2, x, d) || !sp_spp (7, x, d) || !sp_spp (61, x, d))
75 : return 0;
76 : }
77 : else
78 : {
79 : ASSERT (SP_NUMB_BITS <= 64);
80 : /* 64-bit primality test
81 : * follows from results by Jaeschke, "On strong pseudoprimes to several
82 : * bases" Math. Comp. 61 (1993) p916 */
83 :
84 453560 : if (!sp_spp (2, x, d) || !sp_spp (3, x, d) || !sp_spp (5, x, d)
85 26236 : || !sp_spp (7, x, d) || !sp_spp (11, x, d) || !sp_spp (13, x, d)
86 26236 : || !sp_spp (17, x, d) || ! sp_spp (19, x, d) || !sp_spp (23, x, d)
87 26236 : || !sp_spp (29, x, d))
88 427324 : return 0;
89 : }
90 :
91 26236 : return 1;
92 : }
93 :
94 : #define CACHE_LINE_SIZE 64
95 :
96 : void *
97 1971866 : sp_aligned_malloc (size_t len)
98 : {
99 : void *ptr, *aligned_ptr;
100 : size_t addr;
101 :
102 1971866 : ptr = malloc (len + CACHE_LINE_SIZE);
103 1971866 : if (ptr == NULL)
104 0 : return NULL;
105 :
106 1971866 : addr = (size_t)ptr;
107 1971866 : addr = CACHE_LINE_SIZE - (addr % CACHE_LINE_SIZE);
108 1971866 : aligned_ptr = (void *)((char *)ptr + addr);
109 :
110 1971866 : *( (void **)aligned_ptr - 1 ) = ptr;
111 1971866 : return aligned_ptr;
112 : }
113 :
114 : void
115 1971866 : sp_aligned_free (void *newptr)
116 : {
117 : void *ptr;
118 :
119 1971866 : if (newptr == NULL)
120 0 : return;
121 1971866 : ptr = *( (void **)newptr - 1 );
122 1971866 : free (ptr);
123 : }
|