My Project
prime.cc
Go to the documentation of this file.
1 /*****************************************
2  * Computer Algebra System SINGULAR *
3  *****************************************/
4 /*
5  * simple prime test for int
6  */
7 
8 #include <cmath>
9 #include "misc/prime.h"
10 #include "factory/cf_primes.h"
11 
12 static int iiIsPrime0(unsigned p) /* brute force !!!! */
13 {
14  unsigned i,j=0 /*only to avoid compiler warnings*/;
15  if (p<=32749) // max. small prime in factory
16  {
17  int a=0;
18  int e=cf_getNumSmallPrimes()-1;
19  i=e/2;
20  do
21  {
23  if (p==j) return p;
24  if (p<j) e=i-1;
25  else a=i+1;
26  i=a+(e-a)/2;
27  } while ( a<= e);
28  if (p>j) return j;
29  else return cf_getSmallPrime(i-1);
30  }
31  unsigned end_i=cf_getNumSmallPrimes()-1;
32  unsigned end_p=(unsigned)sqrt((double)p);
33 restart:
34  for (i=0; i<end_i; i++)
35  {
37  if ((p%j) == 0)
38  {
39  if (p<=32751) return iiIsPrime0(p-2);
40  p-=2;
41  goto restart;
42  }
43  if (j > end_p) return p;
44  }
45  if (i>=end_i)
46  {
47  while(j<=end_p)
48  {
49  j+=2;
50  if ((p%j) == 0)
51  {
52  if (p<=32751) return iiIsPrime0(p-2);
53  p-=2;
54  goto restart;
55  }
56  }
57  }
58  return p;
59 }
60 
61 int IsPrime(int p) /* brute force !!!! */
62 {
63  if (p == 0) return 0;
64  else if (p == 1) return 1/*1*/;
65  else if ((p == 2)||(p==3)) return p;
66  else if (p < 0) return 2; //(iiIsPrime0((unsigned)(-p)));
67  else if ((p & 1)==0) return iiIsPrime0((unsigned)(p-1));
68  return iiIsPrime0((unsigned)(p));
69 }
70 
int i
Definition: cfEzgcd.cc:132
int p
Definition: cfModGcd.cc:4080
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
access to prime tables
int j
Definition: facHensel.cc:110
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
int IsPrime(int p)
Definition: prime.cc:61
static int iiIsPrime0(unsigned p)
Definition: prime.cc:12