File: sp.c

package info (click to toggle)
gmp-ecm 7.0.4+ds-5
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 4,728 kB
  • sloc: asm: 36,431; ansic: 34,057; xml: 885; python: 799; sh: 698; makefile: 348
file content (123 lines) | stat: -rw-r--r-- 2,999 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
/* sp.c - "small prime" functions that don't need to be inlined

Copyright 2005, 2006, 2007, 2008, 2009, 2010 Dave Newman, Jason Papadopoulos,
Alexander Kruppa, Paul Zimmermann.

The SP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The SP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the SP Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
MA 02110-1301, USA. */

#include <stdio.h> /* for stderr */
#include <stdlib.h>
#include "sp.h"

/* Test if m is a base "a" strong probable prime */

static int
sp_spp (sp_t a, sp_t m, sp_t d)
{
  sp_t r, s, t, e;

  /* since this routine is called only with m >= SP_MIN >= 2^30,
     and a <= 61, we cannot have m = a */
  ASSERT(a < m);
	
  /* Set e * 2^s = m-1, e odd */
  for (s = 0, e = m - 1; !(e & 1); s++, e >>= 1);

  t = sp_pow (a, e, m, d);
	
  if (t == 1)
    return 1;
	
  for (r = 0; r < s; r++)
    {
      if (t == m - 1)
        return 1;

      t = sp_sqr (t, m, d);
    }
	
  return 0;
}

/* Test if x is a prime, return 1 if it is. Note this only works on sp's, 
   i.e. we need the top bit of x set.
   Assume x is odd and x >= SP_MIN.
 */
int
sp_prime (sp_t x)
{
  sp_t d;

  ASSERT (x & 1);
  ASSERT (x >= SP_MIN);

  sp_reciprocal (d, x);
  
  if (SP_NUMB_BITS <= 32)
    {
      /* 32-bit primality test
       * See http://primes.utm.edu/prove/prove2_3.html */
      
      if (!sp_spp (2, x, d) || !sp_spp (7, x, d) || !sp_spp (61, x, d))
        return 0;
    }
  else
    {
      ASSERT (SP_NUMB_BITS <= 64);
      /* 64-bit primality test
       * follows from results by Jaeschke, "On strong pseudoprimes to several
       * bases" Math. Comp. 61 (1993) p916 */
      
      if (!sp_spp (2, x, d) || !sp_spp (3, x, d) || !sp_spp (5, x, d)
        || !sp_spp (7, x, d) || !sp_spp (11, x, d) || !sp_spp (13, x, d)
        || !sp_spp (17, x, d) || ! sp_spp (19, x, d) || !sp_spp (23, x, d)
        || !sp_spp (29, x, d))
          return 0;
    }
  
  return 1;
}

#define CACHE_LINE_SIZE 64

void *
sp_aligned_malloc (size_t len)
{
  void *ptr, *aligned_ptr;
  size_t addr;

  ptr = malloc (len + CACHE_LINE_SIZE);
  if (ptr == NULL)
    return NULL;

  addr = (size_t)ptr;				
  addr = CACHE_LINE_SIZE - (addr % CACHE_LINE_SIZE);
  aligned_ptr = (void *)((char *)ptr + addr);

  *( (void **)aligned_ptr - 1 ) = ptr;
  return aligned_ptr;
}

void
sp_aligned_free (void *newptr) 
{
  void *ptr;

  if (newptr == NULL) 
    return;
  ptr = *( (void **)newptr - 1 );
  free (ptr);
}