Быстрая проверка натурального числа на простоту

int Prime(unsigned long a)
{
   unsigned long i;
   if (a == 2)
      return 1;
   if (a == 0 || a == 1 || a % 2 == 0)
      return 0;
   for(i = 3; i*i <= a && a % i; i += 2)
      ;
   return i*i > a;
}

int Prime(unsigned long a)
{
   unsigned long i, j, bound;
   if (a == 0 || a == 1)
      return 0;
   if (a == 2 || a == 3 || a == 5)
      return 1;
   if (a%2 == 0 || a%3 == 0 || a%5 == 0)
      return 0;
   bound = sqrt((double)a);
   i = 7; j = 11;
   while (j <= bound && a%i && a%j)
   {
       i += 6; j += 6;
   }
   if (j <= bound || i <= bound && a%i == 0)
      return 0;
   return 1;
}

int Prime(unsigned long a)
{
   unsigned long i1, i2, i3, i4, i5, i6, i7, i8, bound;
   if (a == 0 || a == 1)
      return 0;
   if (a == 2 || a == 3 || a == 5 || a == 7 || a == 11 || a == 13 || a == 17 || a == 19 || a == 23 || a == 29)
      return 1;
   if (a%2 == 0 || a%3 == 0 || a%5 == 0 || a%7 == 0 || a%11 == 0 || a%13 == 0 || a%17 == 0 || a%19 == 0 || a%23 == 0 || a%29 == 0)
      return 0;
   bound = sqrt((double)a);
   i1 = 31; i2 = 37; i3 = 41; i4 = 43; i5 = 47; i6 = 49; i7 = 53; i8 = 59;
   while (i8 <= bound && a%i1 && a%i2 && a%i3 && a%i4 && a%i5 && a%i6 && a%i7 && a%i8)
   {
       i1 += 30; i2 += 30; i3 += 30; i4 += 30; i5 += 30; i6 += 30; i7 += 30; i8 += 30;
   }
   if (i8 <= bound ||
      i1 <= bound && a % i1 == 0 ||
      i2 <= bound && a % i2 == 0 ||
      i3 <= bound && a % i3 == 0 ||
      i4 <= bound && a % i4 == 0 ||
      i5 <= bound && a % i5 == 0 ||
      i6 <= bound && a % i6 == 0 ||
      i7 <= bound && a % i7 == 0)
         return 0;
   return 1;
}

void Eratosfen(char *a, long n)
{
    long i, j, bound = sqrt((double)n);
    a[0] = a[1] = 0; a[2] = 1;
    for(i = 4; i < n; i += 2)
        a[i] = 0;
    for(i = 3; i < n; i += 2)
        a[i] = 1;
    i = 3;
    while (i <= bound)
    {
        for(j = i << 1; j < n; j += i)
            a[j] = 0;
        i += 2;
        while (i <= bound && !a[i])
            i += 2;
    }
}

void Eratosfen(char *a, long n)
{
    long i, j, d, bound = sqrt((double)n);
    a[0] = a[1] = 0; a[2] = 1;
    for(i = 4; i < n; i += 2)
        a[i] = 0;
    for(i = 3; i < n; i += 2)
        a[i] = 1;
    i = 3;
    while (i <= bound)
    {
        d = i << 1;
        for(j = 3 * i; j < n; j += d)
            a[j] = 0;
        i += 2;
        while (i <= bound && !a[i])
            i += 2;
    }
}

void Eratosfen(char *a, long n)
{
    long i, j, d, bound = sqrt((double)n);
    a[0] = a[1] = 0; a[2] = 1;
    for(i = 4; i < n; i += 2)
        a[i] = 0;
    for(i = 3; i < n; i += 2)
        a[i] = 1;
    i = 3;
    while (i <= bound)
    {
        d = i << 1;
        for(j = i * i; j < n; j += d)
            a[j] = 0;
        i += 2;
        while (i <= bound && !a[i])
            i += 2;
    }
}

__global__ void is_simple(ulong *arr, ulong *info, int idx, ulong sqr, ulong shift, ulong count)
{
    int thIdx = blockIdx.x * blockDim.x + threadIdx.x;
    ulong x;
    
    if(thIdx % 2 == 0)
        x = 3 * thIdx + 5;
    else
        x = 3 * (thIdx - 1) + 7;
    
    ulong temp = arr[idx];
    ulong rez = 1;
    
    for(ulong i = 0; i < count; ++i)
    {       
        rez = rez xor ((temp % x) == 0);
        x += shift; 
    }
    info[thIdx] = rez;
}
 
void is_simple_arr(ulong *in_arr,ulong *out_arr, int dim_arr, float *time)
{
    dim3 blocks(60000);
    dim3 threads(256);
    int coef = 3;
    ulong *devArr,*devInfo;
    cudaEvent_t start,stop;
    
    cudaEventCreate(&start);
    cudaEventCreate(&stop); 
    
    cudaMalloc((void**)&devArr, dim_arr * sizeof(ulong));
    cudaMalloc((void**)&devInfo, blocks.x * threads.x * sizeof(ulong));
    
    device_ptr<ulong> for_reduction(devInfo);
    cudaMemcpy(devArr, in_arr, dim_arr * sizeof(ulong), cudaMemcpyHostToDevice);
    ulong red_rez;
    
    cudaEventRecord(start, 0); 
    for(int i = 0; i < dim_arr; ++i)
    {
        ulong sq = sqrt(in_arr[i]) + 1;
        if(in_arr[i] % 2 != 0 && in_arr[i] % 3 != 0)
        {
            is_simple<<< blocks, threads >>>(devArr, devInfo, i, sq, blocks.x * threads.x * coef, sq/(blocks.x * threads.x * coef) + 1);
            red_rez = reduce(for_reduction, for_reduction + blocks.x * threads.x,999,thrust::minimum<ulong>());
            out_arr[i] = red_rez;           
        }
        else
        {
            out_arr[i] = 0;
        }
    }
    
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(time,start, stop);
}

if(thIdx % 2 == 0)

if(!(thIdx & 1))

if(in_arr[i] % 2 != 0)

if(in_arr[i] & 1)

ulong sq = sqrt(in_arr[i]) + 1;

ulong sq = sqrt(in_arr[i]);

for( ; dig != 0; )
{
     dig >>= 1;
     ++sum;
}