From b1cdfa01d9c283b3d59441a8843137986cb76a67 Mon Sep 17 00:00:00 2001 From: ketrptr Date: Thu, 5 Dec 2024 13:29:17 +0100 Subject: [PATCH] new file: primeCalc/.gitignore new file: primeCalc/Makefile new file: primeCalc/README.md new file: primeCalc/src/prime_ex.c new file: primeCalc/src/sieve.c new file: primeCalc/src/sieve.h --- primeCalc/.gitignore | 3 ++ primeCalc/Makefile | 18 +++++++ primeCalc/README.md | 10 ++++ primeCalc/src/prime_ex.c | 22 ++++++++ primeCalc/src/sieve.c | 109 +++++++++++++++++++++++++++++++++++++++ primeCalc/src/sieve.h | 28 ++++++++++ 6 files changed, 190 insertions(+) create mode 100644 primeCalc/.gitignore create mode 100644 primeCalc/Makefile create mode 100644 primeCalc/README.md create mode 100644 primeCalc/src/prime_ex.c create mode 100644 primeCalc/src/sieve.c create mode 100644 primeCalc/src/sieve.h diff --git a/primeCalc/.gitignore b/primeCalc/.gitignore new file mode 100644 index 0000000..afc5370 --- /dev/null +++ b/primeCalc/.gitignore @@ -0,0 +1,3 @@ +prime.txt +fastOut +debugOut diff --git a/primeCalc/Makefile b/primeCalc/Makefile new file mode 100644 index 0000000..6438eb8 --- /dev/null +++ b/primeCalc/Makefile @@ -0,0 +1,18 @@ +CFLAGSD=-g -ggdb -Wall -Wextra -Werror +CFLAGSF=-Wall -Wextra -Ofast -march=native -mtune=native -funroll-all-loops +FILES=src/sieve.c src/sieve.h src/prime_ex.c + +all: fast + +debug: $(FILES) + gcc $(CFLAGSD) $(FILES) -o debugOut + +fast: $(FILES) + gcc $(CFLAGSF) $(FILES) -o fastOut + +clean: + rm -rf *.o debugOut fastOut prime.txt + + + + diff --git a/primeCalc/README.md b/primeCalc/README.md new file mode 100644 index 0000000..8e22aae --- /dev/null +++ b/primeCalc/README.md @@ -0,0 +1,10 @@ +# Prime Sieve + +Prime sieve using bitfields + +'make' for fast benchmark output + +'make debug' for debugOut + +If program does not finish in appropriate time adjust n_per_clock in src/sieve.c + diff --git a/primeCalc/src/prime_ex.c b/primeCalc/src/prime_ex.c new file mode 100644 index 0000000..457f7dd --- /dev/null +++ b/primeCalc/src/prime_ex.c @@ -0,0 +1,22 @@ +#include "sieve.h" +#include +#include + +int main(int argc, char **argv) { + FILE *ptr_file; + ptr_file = fopen("prime.txt", "w"); + + clock_t t_0 = clock(); + //// prime algorithm START //// + sieve(argc, argv, ptr_file); + + //// prime algorithm END //// + clock_t t_E = clock(); + + fclose(ptr_file); + + float t_delta = (float)(t_E - t_0) / CLOCKS_PER_SEC; + printf("\nTime in Seconds: %.3f\n", t_delta); + + return 0; +} diff --git a/primeCalc/src/sieve.c b/primeCalc/src/sieve.c new file mode 100644 index 0000000..579bb58 --- /dev/null +++ b/primeCalc/src/sieve.c @@ -0,0 +1,109 @@ +/*Author: Simon Schurti + *License: MIT + * + * + * Prime sieve using bitfields & masks + */ + +#include "sieve.h" + +// primes up to n per clock +// may need adjusting on different platform to finish in appropiate time +const WORD n_per_clock = 130; + +#include +#include +#include + +WORD upperLimit = 4e9; +// WORD upperLimit = 542; + +// fast square root approximation function for unsigned ints +// uses Herons Method +// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method +WORD root(WORD n) { + + WORD x = 0x1 << (MSB(n) - 2); // initial estimate 2 orders smaller than n + WORD prev = 0; + while (x != prev) { + prev = x; + x = (x + (n / x)) / 2; + } + return x; +} + +// find next set bit > prime +// this will be next prime as all numbers up to next prime should be eliminated +WORD findNextPrime(WORD *sieve, WORD prime) { + WORD next = prime + 2; + /*printf("word: %lb\n", sieve[WordIdx(next)]);*/ + while (IsSet(sieve[WordIdx(next)], next) == 0) { + // printf("%lu is not prime\n", next); + next += 2; + if (next >= upperLimit) + return 0; + } + return next; +} + +int sieve(int argc, char **argv, FILE *fp) { + int targetTime = 15; + if (argc > 1) { + targetTime = atoi(argv[1]); + } + upperLimit = (targetTime * CLOCKS_PER_SEC) * n_per_clock; + + // upperLimit = 1e6; + + printf("calculating0 up to : %lu\n", upperLimit); + + // 1 bit per prime candidate / 2 beacause all even can be ignored + WORD buffSize = upperLimit / 2 / BitsPerWord; + + // main sieve buffer + WORD *sieve = (WORD *)calloc(buffSize, BitsPerWord); + if (sieve == NULL) { + perror("calloc: "); + exit(1); + } + + // initialize all as prime + + WORD i = buffSize + 1; + while (i--) { + sieve[i] |= ones; + } + + WORD idx; + WORD nextPrime = 3; // 2 can be assumed 3=2+1 + WORD MaxPrime = nextPrime; + WORD primesFound = 1; // 2 is prime + + //// findNextPrime returns 0 if none found in limit + while (nextPrime > 0 && nextPrime < upperLimit) { + primesFound++; + WORD maxFactor = upperLimit / nextPrime; + for (WORD factor = nextPrime; factor <= maxFactor; factor += 2) { + idx = factor * nextPrime; + if (idx >= upperLimit) + break; + /*printf("eliminating %lu because %lu * %lu = %lu\n", idx, factor,*/ + /* nextPrime, factor * nextPrime);*/ + // printf("with mask: %llb\n", BitIdxMask(idx)); + + sieve[WordIdx(idx)] &= + sieve[WordIdx(idx)] ^ + (BitIdxMask(idx)); // set multiple of next prime to 0 + + // printf("new bu[%lu]: %lb\n", WordIdx(idx), sieve[WordIdx(idx)]); + } + // printf("-----> %lu %llb\n", nextPrime, BitIdxMask(nextPrime)); + MaxPrime = nextPrime; + fprintf(fp, "%lu\n", nextPrime); // write prime number to textfile + nextPrime = findNextPrime(sieve, nextPrime); + } + printf("found %lu primes\nmax prime: %lu\n", primesFound, MaxPrime); + free(sieve); + + return 0; +} diff --git a/primeCalc/src/sieve.h b/primeCalc/src/sieve.h new file mode 100644 index 0000000..88a45ac --- /dev/null +++ b/primeCalc/src/sieve.h @@ -0,0 +1,28 @@ +#ifndef SIEVE +#define SIEVE + +#include +#include +#include + +#define WORD uint64_t +#define BYTE uint8_t +#define BitsPerWord 64 +#define ones 0xFFFFFFFFFFFFFFFFuLL + +// find index of word that contains nth bit in total buffer /2 because no evens +#define WordIdx(n) (((n) / 2) / BitsPerWord) +// bit offset in word +#define BitIdx(n) ((((n) / 2) % BitsPerWord)) + +// mask(3) = 0b10 +#define BitIdxMask(n) 0x1ull << (BitIdx(n)) + +#define IsSet(w, b) (w & BitIdxMask(b)) + +// find most significant bit +//__builtin_clz should work with gcc on all platforms +#define MSB(n) (BitsPerWord - 1 - __builtin_clz(n)) +int sieve(int argc, char **argv, FILE *fp); + +#endif // !SIEVE