TheAlgorithms-C/project_euler/Problem 23/sol2.c

154 lines
3.5 KiB
C

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#ifdef _OPENMP
#include <omp.h>
#endif
/**
* Optimization 1 - compute & store abundant numbers once
* into a look-up array
**/
unsigned long MAX_N = 28123;
char *abundant_flags = NULL;
/**
* Returns:
* -1 if N is deficient
* 1 if N is abundant
* 0 if N is perfect
**/
char get_perfect_number(unsigned long N)
{
unsigned long sum = 1;
char ret = 0;
for (unsigned long i = 2; i * i <= N; i++)
{
if (N % i == 0)
{
sum += i;
unsigned long tmp = N / i;
if (tmp != i)
sum += tmp;
}
}
ret = sum == N ? 0 : (sum > N ? 1 : -1);
#ifdef DEBUG
printf("%5lu: %5lu : %d\n", N, sum, ret);
#endif
return ret;
}
/**
* Find the next abundant number after N and not including N
**/
char is_abundant(unsigned long N)
{
// return abundant_flags[N >> 3] & (1 << N % 8) ? 1 : 0;
return abundant_flags[N >> 3] & (1 << (N & 7)) ? 1 : 0;
}
/**
* Find the next abundant number after N and not including N
**/
unsigned long get_next_abundant(unsigned long N)
{
unsigned long i;
for (i = N + 1; !is_abundant(i); ++i)
;
return i;
}
/**
* check if a given number can be represented as a sum
* of two abundant numbers.
* 1 - if yes
* 0 - if not
**/
char is_sum_of_abundant(unsigned long N)
{
for (unsigned long i = get_next_abundant(1); i <= (N >> 1); i = get_next_abundant(i))
if (is_abundant(N - i))
{
#ifdef DEBUG
printf("\t%4lu + %4lu = %4lu\n", i, N - i, N);
#endif
return 1;
}
return 0;
}
int main(int argc, char **argv)
{
unsigned long sum = 0;
if (argc == 2)
MAX_N = strtoul(argv[1], NULL, 10);
/** byte array to store flags to identify abundant numbers
* the flags are identified by bits
**/
abundant_flags = (char *)calloc(MAX_N >> 3, 1);
#ifdef _OPENMP
printf("Using OpenMP parallleization with %d threads\n", omp_get_max_threads());
#else
printf("Not using parallleization!\n");
#endif
clock_t start_time = clock();
/* Loop to set abundant flags */
#ifdef _OPENMP
#pragma omp for schedule(runtime)
#endif
for (unsigned long N = 1; N <= MAX_N; N++)
{
char ret = get_perfect_number(N);
if (ret == 1)
{
// int byte_offset = N % 8, index = N >> 3;
int byte_offset = N & 7, index = N >> 3;
#ifdef _OPENMP
#pragma omp critical
#endif
abundant_flags[index] |= ret << byte_offset;
}
// if (i % 100 == 0)
// printf("... %5lu: %8lu\r", i, sum);
}
clock_t end_time = clock();
double t1 = 1e3 * (end_time - start_time) / CLOCKS_PER_SEC;
printf("Time taken to get abundant numbers: %.4g ms\n", t1);
start_time = clock();
#ifdef _OPENMP
#pragma omp parallel for reduction(+ \
: sum) schedule(runtime)
#endif
for (unsigned long i = 1; i < MAX_N; i++)
{
if (!is_sum_of_abundant(i))
sum += i;
// if (i % 100 == 0)
// printf("... %5lu\n", i);
}
end_time = clock();
#ifdef DEBUG
putchar('\n');
#endif
double t2 = 1e3 * (end_time - start_time) / CLOCKS_PER_SEC;
printf("Time taken for final sum: %.4g ms\nTotal Time taken: %.4g ms\n", t2, t1 + t2);
printf("Sum of numbers that cannot be represented as sum of two abundant numbers : %lu\n", sum);
free(abundant_flags);
// free(abundant_sum_flags);
return 0;
}