mirror of https://github.com/TheAlgorithms/C
148 lines
3.1 KiB
C
148 lines
3.1 KiB
C
/**
|
|
* \file
|
|
* \brief [Problem 401](https://projecteuler.net/problem=401) solution -
|
|
* Sum of squares of divisors
|
|
* \author [Krishna Vedala](https://github.com/kvedala)
|
|
*/
|
|
#include <stdint.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <time.h>
|
|
#define __STDC_FORMAT_MACROS
|
|
#include <inttypes.h>
|
|
#ifdef _OPENMP
|
|
#include <omp.h>
|
|
#endif
|
|
|
|
#define MOD (uint64_t)1e9 /**< modulo limit */
|
|
#define MAX_L 5000 /**< chunk size of array allocation */
|
|
|
|
/**
|
|
* Check if a number is present in given array
|
|
* \param[in] N number to check
|
|
* \param[in] D array to check
|
|
* \param[in] L length of array
|
|
* \returns 1 if present
|
|
* \returns 0 if absent
|
|
*/
|
|
char is_in(uint64_t N, uint64_t *D, uint64_t L)
|
|
{
|
|
uint64_t i;
|
|
for (i = 0; i < L; i++)
|
|
if (D[i] == N)
|
|
return 1;
|
|
return 0;
|
|
}
|
|
|
|
/**
|
|
* Get all integer divisors of a number
|
|
* \param[in] N number to find divisors for
|
|
* \param[out] D array to store divisors in
|
|
* \returns number of divisors found
|
|
*/
|
|
uint64_t get_divisors(uint64_t N, uint64_t *D)
|
|
{
|
|
uint64_t q, r;
|
|
int64_t i, num = 0;
|
|
|
|
if (N == 1)
|
|
{
|
|
D[0] = 1;
|
|
return 1;
|
|
}
|
|
|
|
// search till sqrt(N)
|
|
// because after this, the pair of divisors will repeat themselves
|
|
for (i = 1; i * i <= N + 1; i++)
|
|
{
|
|
r = N % i; // get reminder
|
|
|
|
// reminder = 0 if 'i' is a divisor of 'N'
|
|
if (r == 0)
|
|
{
|
|
q = N / i;
|
|
if (!is_in(i, D, num)) // if divisor was already stored
|
|
{
|
|
D[num] = i;
|
|
num++;
|
|
}
|
|
if (!is_in(q, D, num)) // if divisor was already stored
|
|
{
|
|
D[num] = q;
|
|
num++;
|
|
}
|
|
}
|
|
|
|
if (num == MAX_L) // limit of array reached, allocate more space
|
|
D = (uint64_t *)realloc(D, MAX_L * sizeof(uint64_t) << 1);
|
|
}
|
|
return num;
|
|
}
|
|
|
|
/**
|
|
* compute sum of squares of all integer factors of a number
|
|
* \param[in] N
|
|
* \returns sum of squares
|
|
*/
|
|
uint64_t sigma2(uint64_t N)
|
|
{
|
|
uint64_t sum = 0, L;
|
|
int64_t i;
|
|
uint64_t *D = (uint64_t *)malloc(MAX_L * sizeof(uint64_t));
|
|
|
|
L = get_divisors(N, D);
|
|
for (i = 1; i < L; i++)
|
|
{
|
|
uint64_t DD = (D[i] * D[i]) % MOD;
|
|
sum += DD;
|
|
}
|
|
|
|
free(D);
|
|
return sum % MOD;
|
|
}
|
|
|
|
/**
|
|
* sum of squares of factors of numbers
|
|
* from 1 thru N
|
|
*/
|
|
uint64_t sigma(uint64_t N)
|
|
{
|
|
uint64_t s, sum = 0;
|
|
int64_t i;
|
|
|
|
#ifdef _OPENMP
|
|
// parallelize on threads
|
|
#pragma omp parallel for reduction(+ : sum)
|
|
#endif
|
|
for (i = 0; i <= N; i++)
|
|
{
|
|
s = sigma2(i);
|
|
sum += s;
|
|
}
|
|
return sum % MOD;
|
|
}
|
|
|
|
/** Main function */
|
|
int main(int argc, char **argv)
|
|
{
|
|
uint64_t N = 1000;
|
|
|
|
if (argc == 2)
|
|
N = strtoll(argv[1], NULL, 10);
|
|
else if (N > 2)
|
|
{
|
|
fprintf(stderr, "Wrong number of input arguments!\n");
|
|
printf("Usage:\t ./sol1.c [N=1000]");
|
|
return -1;
|
|
}
|
|
|
|
clock_t start_time = clock();
|
|
uint64_t result = sigma(N);
|
|
double dtime = clock() - start_time;
|
|
|
|
printf("N = %" PRIu64 "\nSum: %" PRIu64 "\n", N, result);
|
|
printf("Time taken: %.4gms\n", dtime * 1e3 / CLOCKS_PER_SEC);
|
|
|
|
return 0;
|
|
}
|