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             : }