Build into dist/, list primes to 1,000,000 with timing, add wiki
This commit is contained in:
@@ -0,0 +1 @@
|
|||||||
|
dist/
|
||||||
@@ -1,36 +1,39 @@
|
|||||||
# Builds prime_avx2 on Linux (gcc) and Windows (MinGW-w64 / MSYS2 gcc).
|
# Builds prime_avx2 into dist/ on Linux (gcc) and Windows (MinGW-w64 / MSYS2 gcc).
|
||||||
# make build
|
# make build
|
||||||
# make run build and run
|
# make run build and run
|
||||||
# make clean remove build output
|
# make clean remove dist/
|
||||||
|
|
||||||
CC ?= gcc
|
CC ?= gcc
|
||||||
CFLAGS ?= -O3 -mavx2 -Wall -Wextra
|
CFLAGS ?= -O3 -mavx2 -Wall -Wextra
|
||||||
LDLIBS += -lm
|
LDLIBS += -lm
|
||||||
|
|
||||||
|
DIST := dist
|
||||||
|
|
||||||
ifeq ($(OS),Windows_NT)
|
ifeq ($(OS),Windows_NT)
|
||||||
EXE := .exe
|
EXE := .exe
|
||||||
RM := del /Q /F
|
|
||||||
# make on Windows may use cmd.exe, which needs backslashes
|
# make on Windows may use cmd.exe, which needs backslashes
|
||||||
RUN := .\prime_avx2$(EXE)
|
RUN := $(DIST)\prime_avx2$(EXE)
|
||||||
RMTARGET := prime_avx2$(EXE)
|
RMDIR := if exist $(DIST) rmdir /S /Q $(DIST)
|
||||||
else
|
else
|
||||||
EXE :=
|
EXE :=
|
||||||
RM := rm -f
|
RUN := ./$(DIST)/prime_avx2
|
||||||
RUN := ./prime_avx2
|
RMDIR := rm -rf $(DIST)
|
||||||
RMTARGET := prime_avx2
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
TARGET := prime_avx2$(EXE)
|
TARGET := $(DIST)/prime_avx2$(EXE)
|
||||||
|
|
||||||
.PHONY: all run clean
|
.PHONY: all run clean
|
||||||
|
|
||||||
all: $(TARGET)
|
all: $(TARGET)
|
||||||
|
|
||||||
$(TARGET): prime_avx2.c
|
$(TARGET): src/prime_avx2.c | $(DIST)
|
||||||
$(CC) $(CFLAGS) $< -o $@ $(LDLIBS)
|
$(CC) $(CFLAGS) $< -o $@ $(LDLIBS)
|
||||||
|
|
||||||
|
$(DIST):
|
||||||
|
mkdir $(DIST)
|
||||||
|
|
||||||
run: $(TARGET)
|
run: $(TARGET)
|
||||||
$(RUN)
|
$(RUN)
|
||||||
|
|
||||||
clean:
|
clean:
|
||||||
-$(RM) $(RMTARGET)
|
-$(RMDIR)
|
||||||
|
|||||||
@@ -0,0 +1,24 @@
|
|||||||
|
# Build and run
|
||||||
|
|
||||||
|
**English**
|
||||||
|
|
||||||
|
The `Makefile` is in the repository root and puts the executable in `dist/`.
|
||||||
|
It needs `gcc` (Linux, or MinGW-w64 / MSYS2 on Windows) and a CPU with AVX2.
|
||||||
|
|
||||||
|
```bash
|
||||||
|
make
|
||||||
|
```
|
||||||
|
|
||||||
|
```bash
|
||||||
|
make run
|
||||||
|
```
|
||||||
|
|
||||||
|
```bash
|
||||||
|
make clean
|
||||||
|
```
|
||||||
|
|
||||||
|
To save the primes and still see the time, redirect only standard output:
|
||||||
|
|
||||||
|
```bash
|
||||||
|
./dist/prime_avx2 > primes.txt
|
||||||
|
```
|
||||||
@@ -0,0 +1,34 @@
|
|||||||
|
# Description
|
||||||
|
|
||||||
|
**English**
|
||||||
|
|
||||||
|
## Purpose
|
||||||
|
|
||||||
|
The program shows how fast a CPU is. How fast it looks depends on how the
|
||||||
|
program is written, so the number is a measure of the CPU and the code
|
||||||
|
together, not of the CPU alone.
|
||||||
|
|
||||||
|
## What it does
|
||||||
|
|
||||||
|
It tests every number from 0 to 1,000,000 with `is_prime_avx2`, which uses
|
||||||
|
AVX2 trial division by 6k ± 1 (four divisors per step). It prints each prime
|
||||||
|
on its own line and writes a summary line to standard error:
|
||||||
|
|
||||||
|
```text
|
||||||
|
Found 78498 primes from 0 to 1000000 in 13.806 ms
|
||||||
|
```
|
||||||
|
|
||||||
|
78,498 is the known count of primes up to 1,000,000, so the summary line also
|
||||||
|
works as a correctness check.
|
||||||
|
|
||||||
|
## Why the time depends on the code
|
||||||
|
|
||||||
|
The same CPU gives very different times for different programs:
|
||||||
|
|
||||||
|
- the algorithm (trial division, a sieve, a segmented sieve),
|
||||||
|
- the instruction set used (scalar code, AVX2, AVX-512),
|
||||||
|
- compiler and flags (`-O3 -mavx2`, MSVC `/O2 /arch:AVX2`),
|
||||||
|
- what is measured (only the search, or also printing to a console).
|
||||||
|
|
||||||
|
Here only the search is timed. Output is collected first and printed
|
||||||
|
afterwards, so console speed does not change the number.
|
||||||
@@ -0,0 +1,10 @@
|
|||||||
|
# prime-numbers-c
|
||||||
|
|
||||||
|
**English**
|
||||||
|
|
||||||
|
A small C program that lists every prime number from 0 to 1,000,000 and
|
||||||
|
reports how long it took.
|
||||||
|
|
||||||
|
- [Description](Description): what the program is for
|
||||||
|
- [Results](Results): measured times on real hardware
|
||||||
|
- [Build and run](Build-and-Run): `make`, `make run`
|
||||||
@@ -0,0 +1,14 @@
|
|||||||
|
# Results
|
||||||
|
|
||||||
|
**English**
|
||||||
|
|
||||||
|
Time to find all primes from 0 to 1,000,000. Only the search is timed.
|
||||||
|
|
||||||
|
| CPU | Result |
|
||||||
|
| --- | --- |
|
||||||
|
| 13th Gen Intel(R) Core(TM) i5-13600K | Found 78498 primes from 0 to 1000000 in 20.977 ms |
|
||||||
|
|
||||||
|
A result is only comparable with another when the code, compiler and flags
|
||||||
|
are the same. To add a result, run the program on your machine (see
|
||||||
|
[Build and run](Build-and-Run)) and add a row with your CPU name and the
|
||||||
|
summary line.
|
||||||
@@ -0,0 +1,4 @@
|
|||||||
|
- [Home](Home)
|
||||||
|
- [Description](Description)
|
||||||
|
- [Results](Results)
|
||||||
|
- [Build and run](Build-and-Run)
|
||||||
@@ -1,59 +1,72 @@
|
|||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <stdbool.h>
|
#include <stdbool.h>
|
||||||
#include <stdint.h>
|
#include <stdint.h>
|
||||||
#include <inttypes.h>
|
#include <inttypes.h>
|
||||||
#include <math.h>
|
#include <math.h>
|
||||||
#include <immintrin.h> // AVX2 intrinsics
|
#include <time.h>
|
||||||
|
#include <immintrin.h> // AVX2 intrinsics
|
||||||
// Returns true if n is prime, using AVX2 trial division by 6k +/- 1.
|
|
||||||
// Supports n up to UINT32_MAX.
|
// Returns true if n is prime, using AVX2 trial division by 6k +/- 1.
|
||||||
//
|
// Supports n up to UINT32_MAX.
|
||||||
// Four divisors are tested per iteration in double precision. A double holds
|
//
|
||||||
// every 32-bit integer exactly and IEEE division is correctly rounded, so
|
// Four divisors are tested per iteration in double precision. A double holds
|
||||||
// floor(n / d) is exact here (a 24-bit float quotient is not).
|
// every 32-bit integer exactly and IEEE division is correctly rounded, so
|
||||||
bool is_prime_avx2(uint64_t n) {
|
// floor(n / d) is exact here (a 24-bit float quotient is not).
|
||||||
if (n <= 1) return false;
|
bool is_prime_avx2(uint64_t n) {
|
||||||
if (n <= 3) return true;
|
if (n <= 1) return false;
|
||||||
if (n % 2 == 0 || n % 3 == 0) return false;
|
if (n <= 3) return true;
|
||||||
if (n > UINT32_MAX) return false; // out of supported range
|
if (n % 2 == 0 || n % 3 == 0) return false;
|
||||||
|
if (n > UINT32_MAX) return false; // out of supported range
|
||||||
const double limit = floor(sqrt((double)n));
|
|
||||||
|
const double limit = floor(sqrt((double)n));
|
||||||
__m256d v_n = _mm256_set1_pd((double)n);
|
|
||||||
__m256d v_limit = _mm256_set1_pd(limit);
|
__m256d v_n = _mm256_set1_pd((double)n);
|
||||||
|
__m256d v_limit = _mm256_set1_pd(limit);
|
||||||
// Divisors 5, 7, 11, 13 (6k +/- 1); each iteration advances by 12.
|
|
||||||
__m128i v_div_i = _mm_setr_epi32(5, 7, 11, 13);
|
// Divisors 5, 7, 11, 13 (6k +/- 1); each iteration advances by 12.
|
||||||
const __m128i v_step = _mm_set1_epi32(12);
|
__m128i v_div_i = _mm_setr_epi32(5, 7, 11, 13);
|
||||||
|
const __m128i v_step = _mm_set1_epi32(12);
|
||||||
for (double base = 5.0; base <= limit; base += 12.0) {
|
|
||||||
__m256d v_d = _mm256_cvtepi32_pd(v_div_i);
|
for (double base = 5.0; base <= limit; base += 12.0) {
|
||||||
__m256d v_q = _mm256_floor_pd(_mm256_div_pd(v_n, v_d));
|
__m256d v_d = _mm256_cvtepi32_pd(v_div_i);
|
||||||
__m256d v_rem = _mm256_sub_pd(v_n, _mm256_mul_pd(v_q, v_d));
|
__m256d v_q = _mm256_floor_pd(_mm256_div_pd(v_n, v_d));
|
||||||
|
__m256d v_rem = _mm256_sub_pd(v_n, _mm256_mul_pd(v_q, v_d));
|
||||||
__m256d is_zero = _mm256_cmp_pd(v_rem, _mm256_setzero_pd(), _CMP_EQ_OQ);
|
|
||||||
__m256d in_range = _mm256_cmp_pd(v_d, v_limit, _CMP_LE_OQ);
|
__m256d is_zero = _mm256_cmp_pd(v_rem, _mm256_setzero_pd(), _CMP_EQ_OQ);
|
||||||
|
__m256d in_range = _mm256_cmp_pd(v_d, v_limit, _CMP_LE_OQ);
|
||||||
// A lane only counts if it divides n and the divisor is <= sqrt(n)
|
|
||||||
if (_mm256_movemask_pd(_mm256_and_pd(is_zero, in_range)) != 0)
|
// A lane only counts if it divides n and the divisor is <= sqrt(n)
|
||||||
return false;
|
if (_mm256_movemask_pd(_mm256_and_pd(is_zero, in_range)) != 0)
|
||||||
|
return false;
|
||||||
v_div_i = _mm_add_epi32(v_div_i, v_step);
|
|
||||||
}
|
v_div_i = _mm_add_epi32(v_div_i, v_step);
|
||||||
|
}
|
||||||
return true;
|
|
||||||
}
|
return true;
|
||||||
|
}
|
||||||
int main(void) {
|
|
||||||
// Known large prime: 1,000,000,007
|
#define LIMIT 1000000u
|
||||||
uint64_t num = 1000000007;
|
|
||||||
|
// Seconds elapsed between two timespec_get(TIME_UTC) samples.
|
||||||
printf("Checking if %" PRIu64 " is prime using AVX2...\n", num);
|
static double elapsed(const struct timespec *a, const struct timespec *b) {
|
||||||
|
return (double)(b->tv_sec - a->tv_sec) + (double)(b->tv_nsec - a->tv_nsec) / 1e9;
|
||||||
if (is_prime_avx2(num))
|
}
|
||||||
printf("%" PRIu64 " is a PRIME number.\n", num);
|
|
||||||
else
|
int main(void) {
|
||||||
printf("%" PRIu64 " is a COMPOSITE number.\n", num);
|
static uint32_t primes[LIMIT / 2]; // pi(1e6) = 78498, so this is ample
|
||||||
|
size_t count = 0;
|
||||||
return 0;
|
struct timespec t0, t1;
|
||||||
}
|
|
||||||
|
// Time the computation only; printing is done afterwards.
|
||||||
|
timespec_get(&t0, TIME_UTC);
|
||||||
|
for (uint32_t n = 0; n <= LIMIT; n++)
|
||||||
|
if (is_prime_avx2(n)) primes[count++] = n;
|
||||||
|
timespec_get(&t1, TIME_UTC);
|
||||||
|
|
||||||
|
for (size_t i = 0; i < count; i++)
|
||||||
|
printf("%" PRIu32 "\n", primes[i]);
|
||||||
|
|
||||||
|
fprintf(stderr, "Found %zu primes from 0 to %u in %.3f ms\n",
|
||||||
|
count, LIMIT, elapsed(&t0, &t1) * 1000.0);
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
Reference in New Issue
Block a user