Приглашаю вас в небольшое приключение выходного дня, в котором никто никому ничего не будет доказывать. Мы просто будем реализовывать один и тот же несложный алгоритм, разыскивающий простые числа в некотором диапазоне, на нескольких языках программирования: C, C++, Scheme и Python - и смотреть, что с этим кодом могут сделать современные оптимизирующие компиляторы. В процессе приключения мы увидим, что «динамический» не означает «совсем уж медленный», и посмотрим на приёмы программирования на Scheme, что, как мне кажется, можно сравнить с путешествием на экзотический остров.
Начало
Началось всё одним томным вечером, когда мне на глаза попалось видео, в котором автор сравнивал три языка программирования: Python, C и Ассемблер IA-32 - по продуктивности программиста и по эффективности полученного кода. Задача решается несложная: подсчитать количество простых чисел не больших 250000.
На мой взгляд, такая оценка - комбинация продуктивности и производительности - хорошо подходит для сравнения языков программирования. Поэтому мне стало интересно, какое место среди этих трёх языков займёт Scheme. При этом, мне хотелось проверить, насколько быстро Scheme позволит написать более сложный алгоритм, который позволит получить решение за меньшее время.
Scheme - язык с динамической типизацией, сборкой мусора, динамическими контекстами и замыканиями, от программ на котором не ожидаешь особой прыти. Ограничив перебор потенциальных делителей простыми числами, я рассчитывал получить код, работающий раз в 10 быстрее, чем программа на Python в видео. Вот asciicast реализации этой затеи.
Результат мягко говоря, удивил (уважаемый читатель, не торопись писать комментарий о sqrt(n); для большего драматического эффекта, я оставил эту деталь напоследок).
time scheme --script primes.scm
22044
real 0m1.408sТо есть, Chez Scheme успевает разобрать исходный текст, скомпилировать код (компилятор написан на Scheme) и выполнить его быстрее, чем прокрутит свои циклы полученный из исходника на C код, который даже память не трогает? «Да, ну, не может быть!» - подумал я: «Просто у автора видео медленный процессор.» Естественно, возникает желание проверить всё на одной машине, в данном случае на виртуальной в облаке Яндекс с такой вот конфигурацией.
|
|
|
|
|
|
|
|
|
|
Все исходные тексты, использованные дальше, выложены здесь
Итерация 0: Не верю!
Исходники на C и Python, взяты из ролика. Я полагаю, что все знакомы с Python и C, поэтому их можно не пояснять.
primes-0.c и primes-0.py
Hidden text
// primes-0.c
#include <stdio.h>
int is_prime(int n) {
for(int p = 2; p <= n/2; p++)
if (!(n % p))
return 0;
return 1;
}
void main() {
int n_primes = 0;
for(int n = 2; n < 250001; n++)
n_primes += is_prime(n);
printf("%d\n", n_primes);
}# primes-0.py
def is_prime(n):
for p in range(2, n//2 + 1):
if (not (n%p)):
return 0
return 1
n_primes = 0
for i in range(2, 250001):
n_primes += is_prime(i)
print(str(n_primes))Со Scheme, вероятно, знакомо меньше читателей, поэтому я буду объяснять, как устроен код.
; primes-0.scm
(define (prime? n P) (or (null? P)
(> (car P) (quotient n 2))
(and (positive? (remainder n (car P)))
(prime? n (cdr P)))))
(define primes (cons 2 '()))
(do ((n 3 (+ 2 n))
(P primes (if (not (prime? n primes)) P (begin (set-cdr! P (cons n '())) (cdr P)))))
((> n 250000) (display (length primes)) (newline)))Код primes-0.scm состоит из двух частей.
Часть первая. Процедура prime? возвращает #f (ложь), если находит простой делитель числа n в списке P. В ней использованы списки и логические конструкции Scheme.
Конструкции (or E1 ... En) и (and E1 ... En) работают по сокращённой схеме (short-circuiting). В C бы им соответствовали конструкции E1 || ... || En и E1 && ... && En.
Списки в Lisp со стародавних времён представляют cons-парами. Cons от constructor: в классическом определении символьных выражений (s-выражений), конструирование пары - единственный способ построения неатомарных значений, других конструкторов нет, отсюда и название. Пару конструирует процедура (примитивная, но обычно это не уточняют, в Scheme все процедуры равноправны) cons.
Список элементов - это пара (элемент, список остальных элементов). Из списка L можно получить новый список добавлением элемента e в голову: (cons e L). Список [1, 2, 3, 4] создаётся выражением (cons 1 (cons 2 (cons 3 (cons 4 '())))). Здесь '() - специальное атомарное значение, обозначающее пустой список.
Пары разбирают на первую и вторую компоненты процедуры car и cdr. Названия исторически связаны со способом представления пар в первой реализации Lisp на машине IBM 704. car - contents of the address part of the register, cdr - contents of the decrement part of the register. Под регистром подразумевается ячейка памяти. По одной ячейке на каждую cons-пару - Морфеус увидел бы в этом предназначение.
Если L - список, то (car L) возвращает первый элемент списка L, (cdr L) - остаток списка. Если L = [1, 2, 3, 4], то (car L) = 1 и (cdr L) = [2, 3, 4] = (cons 2 (cons 3 (cons 4 '()))).
Это понадобится чуть позже, но напишу об этом сразу. Scheme (как и многие другие диалекты Lisp) позволяет изменять значения полей cons-ячеек процедурами set-car! и set-cdr!.
Приняв во внимание вышесказанное, код prime?
; primes-0.scm
(define (prime? n P) (or (null? P)
(> (car P) (quotient n 2))
(and (positive? (remainder n (car P)))
(prime? n (cdr P)))))можно прочитать так: (prime? n P) - истина,
если
(null? P)- список закончился;или если не закончился, то если
(> (car P) (quotient n 2))- очередной простое число достаточно большое, и дальше проверять смысла нет;или если смысл проверять есть, то:
если
(positive? (remainder n (car P)))- текущее простое число не делитn,и тогда, если
(prime? n (cdr P))- в остатке списка простых чисел не найдётся делителя дляn.
Стандарт Scheme гарантирует, что в этом случае вызов (prime? n (cdr P)) будет хвостовым, то есть, стек кадров активации процедур не будет расти. Поэтому можно считать, что prime? - реализация цикла.
Вторая часть последовательно перебирает целые числа n от 3 до 250000 c шагом 2 и, обнаружив простое число, добавляет его в конец списка простых чисел primes. Переменная P - курсор, указывающий на этот конец. Реализован этот процесс циклом do
(do ((V1 I1 U1) ... (Vn In Un))
(C E1 ... Ek)
B1
...
Bm)Работает это так:
В области видимости
doсоздаются переменныеV1, ...,Vn. КаждаяVjинициализируется значением соответствующего выраженияIj.На каждой итерации цикла
Vjобновляется значением выраженияUj. Цикл завершается, когда выражениеСстановится истинным, перед выходом из цикла исполняются выраженияE1, ...,Ek. Значением всего выраженияdoстановится результатEk.На каждой итерации выполняются выражения в теле цикла
B1, ...,Bn.
В do в тексте primes-0.scm
(do ((n 3 (+ 2 n))
(P primes (if (not (prime? n primes)) P (begin (set-cdr! P (cons n '())) (cdr P)))))
((> n 250000) (display (length primes)) (newline)))вероятно, самым мудрёным является выражение для обновления переменной P, указывающей на последнюю cons-пару списка primes:
(if (not (prime? n primes)) P (begin (set-cdr! P (cons n '())) (cdr P)))В нём проверяется простота n. Если n не простое, то возвращается значение P без изменений. Если простое, то во второй элемент cons-пары P записывается ссылка на новую cons-пару, содержащую атомарный список [n]. Эта ссылка возвращается последним выражением блока begin - (cdr P) - и становится новым значением переменной P. Инвариант цикла: P указывает на последнюю cons-пару списка primes - сохраняется.
Посмотрим, насколько быстро это всё работает.
$ gcc primes-0.c && time ./a.out
22044
real 0m3.854sАвтор ролика грубо пошутил над питонистами, забыв о том, что они вооружены JIT-компилятором PyPy, который отлично справляется в данном случае.
$ time pypy primes-0.py
22044
real 0m4.859sScheme ...
$ time scheme --script primes-0.scm
22044
real 0m1.468s.. в 2.5 раза быстрее C. «Как же так? Не верю! Что происходит!? В Debian-версии GCC отключены какие-то оптимизации? Надо включить всё.» - вертелось у меня в голове.
$ gcc -O3 primes-0.c && time ./a.out
22044
real 0m2.769s«Да что же это такое!?» - вопрошал мой уже слегка воспалённый разум. Но раз в ход пошла агрессивная оптимизация, то, для полноты картины, нужно применить её и для программы на Scheme.
Во-первых, нужно заменить всю арифметику на арифметику с fixnum-ами. В Lisp fixnum-ами называют целые числа, работа с которыми гарантированно осуществляется напрямую, без boxing-а, то есть, не по ссылкам.
primes-fx-0.scm
Hidden text
; primes-fx-0.scm
(define (prime? n P) (or (null? P)
(fx> (car P) (fxquotient n 2))
(and (fxpositive? (fxremainder n (car P)))
(prime? n (cdr P)))))
(define primes (cons 2 '()))
(do ((n 3 (fx+ 2 n))
(P primes (if (prime? n primes) (begin (set-cdr! P (cons n '())) (cdr P)) P)))
((fx> n 250000) (display (length primes)) (newline)))Во-вторых, нужно разрешить компилятору пуститься во все тяжкие.
$ time scheme --optimize-level 3 --script primes-fx-0.scm
22044
real 0m0.538sОптимизированный код на C, не работающий с памятью, оказался в 5 раз медленней оптимизированного кода на Scheme (языке с динамической типизацией, сборкой мусора), вовсю использующим для своей работы динамическую структуру данных и с циклами в виде рекурсии!?
Вечер стремительно терял свою томность.
Итерация 1: Алгоритмизируй это!
Я был готов к тому, что алгоритмы побеждают грубую силу. Но я не ожидал такого результата (ведь, перебор шёл не до sqrt(n)). Мне захотелось выяснить, насколько компилятор Chez Scheme хорош.
Нужно реализовать похожие алгоритмы на Python, C++ и C, и сравнить.
На Python и C++ это довольно просто сделать. Списки в Python позволяют быстро добавлять элементы в конец. В стандартном C++ есть векторы с более плотной упаковкой значений в памяти и амортизированным добавлением элементов в конец (в этой части моего приключения C++ и появился по этой причине).
Код на C++ и Python должен быть понятен.
python-1.cpp и python-1.py
Hidden text
// primes-1.cpp
#include <iostream>
#include <vector>
using std::vector;
static bool is_prime(const long n, const vector<long> &P) {
for (auto p : P) {
if (p > n / 2) return true;
if (!(n % p)) return false;
}
return true;
}
int main(int argc, char *argv[]) {
if (argc < 2) {
std::cerr << "Specify limit" << std::endl;
exit(1);
}
const unsigned long limit = std::stol(argv[1]);
vector<long> primes = {2};
for(long n = 3; n <= limit; n += 2) {
if (is_prime(n, primes)) {
primes.push_back(n);
}
}
std::cout << primes.size() << std::endl;
}# primes-1.py
import sys
def is_prime(n, P):
for p in P:
if p > n//2: return True
if not (n%p): return False
return True
limit = int(sys.argv[1])
primes = [2]
for n in range(3, limit + 1, 2):
if is_prime(n, primes): primes.append(n)
print(len(primes))Код на Scheme отличается от primes-fx-0.scm тем, что граница поиска простых чисел считывается из командной строки выражением: (string->number (cadr (command-line))). В котором загадочной может показаться процедура cadr. Обычно её определяют так:
(define (cadr L) (car (cdr L)))cadr извлекает первую компоненту из второй cons-пары списка, то есть, второй элемент списка.
primes-1.scm
Hidden text
; primes-1.scm
(define (prime? n P)
(or (null? P)
(fx> (car P) (fxquotient n 2))
(and (fxpositive? (fxremainder n (car P)))
(prime? n (cdr P)))))
(define N (string->number (cadr (command-line))))
(define primes (cons 2 '()))
(do ((n 3 (fx+ 2 n))
(P primes (if (prime? n primes) (begin (set-cdr! P (cons n '())) (cdr P)) P)))
((fx> n N) (display (length primes)) (newline)))Программа на C - калька программы на Scheme с одним отличием. В процессе подсчёта найденных простых чисел в функции counting_free освобождается занятая память. Конкретно в этой программе освобождать память не нужно: ядро ОС освободит ресурсы завершённого процесса. Но в Python и Scheme есть сборщики мусора. В программе на C++ будет вызван деструктор вектора. Поэтому постараемся и на C быть аккуратными.
primes-1.c
Hidden text
// primes-1.c
#include <stdio.h>
#include <stdlib.h>
typedef struct cons {
long p;
struct cons *next;
} Cons;
static Cons *cons(const long v) {
Cons *const c = malloc(sizeof(Cons));
c->p = v;
c->next = NULL;
return c;
}
static int is_prime(const long n, const Cons *const list) {
const Cons *P = list;
for (const Cons *P = list; P; P = P->next) {
if (P->p > n/1) return 1;
if (n % P->p == 0) return 0;
}
return 1;
}
static long counting_free(Cons *const list) {
long n = 0;
const Cons *P = list;
while(P) {
const Cons *const next = P->next;
free((void *) P);
n++;
P = next;
}
return n;
}
int main(int argc, char *argv[]) {
if (argc < 2) {
fprintf(stderr, "Specify limit\n");
return 1;
}
const long N = atol(argv[1]);
Cons *const primes = cons(2);
Cons *P = primes;
for(long n = 3; n <= N; n += 2) {
if (is_prime(n, primes)) {
P->next = cons(n);
P = P->next;
}
}
printf("%ld\n", counting_free(primes));
return 0;
}
Запускаем наших претендентов.
$ time scheme --script primes-1.scm 250000
22044
real 0m1.156s
$ g++ primes-1.cpp && time ./a.out 250000
22044
real 0m1.231s
$ gcc primes-1.c && time ./a.out 250000
22044
real 0m0.972s
$ time pypy primes-1.py 250000
22044
real 0m0.505s«Да вы издеваетесь! PyPy!?» PyPy оказался резвее всех, завершив вычисление почти в 2 раза быстрее ближайшего конкурента. Программы на C++ и С не слишком далеко оторвались от программы на Scheme. Включение 3-го уровня оптимизации сделало программу на C++ в 2.5 раза быстрее (на самом деле, неожиданный и неприятный сюрприз), но значительно превзойти результат PyPy это не позволило. А программа на C оказалась медленней программы на Scheme.
$ gcc -O3 primes-1.c && time ./a.out 250000
22044
real 0m0.862s
$ time scheme --optimize-level 3 --script primes-1.scm 250000
22044
real 0m0.514s
$ g++ -O3 primes-1.cpp && time ./a.out 250000
22044
real 0m0.459s«Да что же это такое!?»
Итерация 2: Дединамизируй это!
Знакомые, посмотрев на безобразие, предположили, что «это» - работа с кучей, и предложили избавиться от malloc-ов. «Что же делать? Где там мой page size?» - подумал я и принялся за дело.
Первый на этой итерации вариант программы на C реализует всё тот же алгоритм, но хранит обнаруженные простые числа в массиве, размер которого увеличивается по необходимости на одну страницу памяти обращением к realloc. Единственная хитрость в коде - подсказка GCC о более вероятном значении выражения A->cursor = A->len при помощи __builtin_expect.
primes-2.c
Hidden text
// primes-2.c
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
typedef struct {
long *restrict buf;
long len;
long cursor;
} Array;
static Array make_array() {
const long pagelen = sysconf(_SC_PAGESIZE) / sizeof(long);
return (Array) {
.buf = malloc(pagelen * sizeof(long)),
.len = pagelen,
.cursor = 0 };
}
#define unlikely(X) __builtin_expect((X), 0)
static void push(const long val, Array *restrict const A) {
if (unlikely(A->cursor >= A->len)) {
A->len += sysconf(_SC_PAGESIZE) / sizeof(long);
A->buf = realloc(A->buf, A->len * sizeof(long));
}
A->buf[A->cursor++] = val;
}
static long counting_free(Array *restrict const A) {
free(A->buf);
A->buf = NULL;
return A->cursor;
}
static int is_prime(const long n, const Array *restrict const P) {
for (int i = 0; i < P->cursor; i++) {
const long p = P->buf[i];
if (p > n / 2) return 1;
if (n % p == 0) return 0;
}
return 1;
}
int main(int argc, char *argv[]) {
if (argc < 2) {
fprintf(stderr, "Specify limit\n");
return 1;
}
const long N = atol(argv[1]);
Array primes = make_array();
push(2, &primes);
for(long n = 3; n <= N; n += 2)
if (is_prime(n, &primes)) push(n, &primes);
printf("%ld\n", counting_free(&primes));
return 0;
}
Такая стратегия работы с памятью позволяет libc использовать механизм виртуальной памяти для релокации массива, чтобы избежать явного копирования данных на новое место. Вместо копирования достаточно перенастраивать таблицу трансляции виртуальных адресов процесса системным вызовом mremap. Что и происходит:
$ gcc primes-2.c && strace ./a.out 250000
...
mremap(0x7f60f2ce0000, 172032, 176128, MREMAP_MAYMOVE) = 0x7f60f2ce0000
mremap(0x7f60f2ce0000, 176128, 180224, MREMAP_MAYMOVE) = 0x7f60f2ce0000
...Снизить интенсивность запросов на выделение памяти можно и в программе на Scheme. Современные Lisp-ы (Chez Scheme не исключение), поддерживают разные виды массивов. Но в стандартной библиотеке Chez нет процедур для эффективного изменения размеров массивов: нужно вручную создавать новый массив с большим размером и копировать в него данные.
Однко можно применить чуть более сложную структуру данных: список fixnum-массивов (fxvector), где каждый вектор размером со страницу. Для каждой такой страницы нужно хранить курсор, указатель на очередной свободный элемент массива. Курсором может быть нулевой элемент вектора. Сказано - сделано.
primes-2.scm
Hidden text
; primes-2.scm
(define page-length 512)
(define (cursor P) (fxvector-ref P 0))
(define (cursor-set! P c) (fxvector-set! P 0 c))
(define (make-page n)
(let ((P (make-fxvector page-length)))
(cursor-set! P 2)
(fxvector-set! P 1 n)
P))
(define (push! n P)
(let ((A (car P))
(c (cursor (car P))))
(if (fx< c (fxvector-length A))
(begin (fxvector-set! A c n)
(cursor-set! A (fx1+ c))
P)
(begin (set-cdr! P (cons (make-page n) '()))
(cdr P)))))
(define (prime? n pages)
(let ((l (fxquotient n 2)))
(let next ((P pages))
(or (null? P)
(let ((A (car P))
(c (cursor (car P))))
(let loop ((i 1))
(if (fx< i c)
(let ((p (fxvector-ref A i)))
(or (fx> p l)
(and (fxpositive? (fxremainder n p))
(loop (fx1+ i)))))
(next (cdr P)))))))))
(define (count P)
(do ((p P (cdr p))
(S 0 (fx+ S (cursor (car p)) -1)))
((null? p) S)))
(let ((N (string->number (cadr (command-line))))
(primes (cons (make-page 2) '())))
(do ((n 3 (fx+ 2 n))
(P primes (if (prime? n primes) (push! n P) P)))
((> n N) (display (count primes)) (newline))))
Первые три процедуры в primes-2.scm задают «архитектуру» страницы.
(define (cursor P) (fxvector-ref P 0))
(define (cursor-set! P c) (fxvector-set! P 0 c))
(define (make-page n)
(let ((P (make-fxvector page-length)))
(cursor-set! P 2)
(fxvector-set! P 1 n)
P))fxvector-ref и fxvector-set! - процедуры для чтения и записи элементов вектора. Процедура make-fxvector создаёт вектор запрошенного размера, инициализируя его нулями.
Конструкция (let ((V1 E1) ... (VN EN)) B1 ... BK) используется для определения переменных в локальной для выражений тела (Bi) области видимости. По смыслу эта конструкция эквивалентна выражению
((lambda (V1 ... VN) B1 ... BK) E1 ... EN)Этот код создаёт анонимную процедуру, параметры которой названы нужным образом, и вызывает её со значениями выражений, которые связываются с соответствующими именам. В резальтате Bi выполняются с нужными значениями переменных Vi.
Процедура push! реализована бесхитростно.
(define (push! n P)
(let ((A (car P))
(c (cursor (car P))))
(if (fx< c (fxvector-length A))
(begin (fxvector-set! A c n)
(cursor-set! A (fx1+ c))
P)
(begin (set-cdr! P (cons (make-page n) '()))
(cdr P)))))Она предполагает, что P указывает на последнюю cons-пару в списке страниц с простыми числами. Если на этой странице осталось место, то на это место записывается число n и курсор сдвигается. Если места нет, создаётся новая страница, в начало которой записывается n, эта страница упаковывается в новую cons-пару, которая «подшивается» в конец списка, и ссылка на этот новый конец списка возвращается.
В процедуре prime? применён другой вариант let: (let N ((V1 E1) ... (VN EN)) B1 ... BK). По смыслу он эквивалентен выражению:
((rec N (lamda (V ...) B ...)) E ...)(rec N E) - оператор неподвижной точки, то есть, он вычисляет выражение E, в окружении, в котором E - значение переменной N.
Например (rec ones (stream-cons 1 ones) будет потоком из единиц:
> (stream->list (stream-take 10 (rec ones (stream-cons 1 ones))))
(1 1 1 1 1 1 1 1 1 1)Или, например, используя rec, можно классическим функциональным способом определить последовательность Фибоначчи
> (stream->list
(stream-take 20
(rec fibs (stream-cons 0 (stream-cons 1 (stream-map + fibs (stream-cdr fibs)))))))
(0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181)Таким образом, выражение
((rec N (lamda (V1 ... VN) B1 ... BK)) E1 ... EN)задаёт процедуру, которая может рекурсивно вызывать себя по имени N, и применяет её к некоторым аргументам. Чувствую, что ещё один пример не помешает. Традиционный пример: факториал.
> ((rec loop (lambda (n r) (if (zero? n) r (loop (- n 1) (* n r))))) 10 1)
3628800Или, эквивалентно, в форме let:
(let loop ((n 10) (r 1)) (if (zero? n) r (loop (- n 1) (* n r))))Вместе с хвостовыми вызовами такой вариант let позволяет описывать циклы со сложной логикой работы, множественными точками передачи управления и выходами.
(define (prime? n pages)
(let ((l (fxquotient n 2)))
(let next ((P pages))
(or (null? P)
(let loop ((i 1))
(if (fx< i (cursor (car P)))
(let ((p (fxvector-ref (car P) i)))
(or (fx> p l)
(and (fxpositive? (fxremainder n p))
(loop (fx1+ i)))))
(next (cdr P))))))))В prime? два цикла: next по списку страниц с простыми числами, и loop по простым числам внутри страницы. Первый выход из next происходит, когда заканчивается список - (null? P) возвращает истину, и эта истина становится результатом внешнего or.
Если же в списке есть страницы, то с первой из них, (car P), начинает работу цикл loop. Если числа на странице закончились - условие (fx< i (cursor (car P)) стало ложным, то происходит переход на следующую итерацию цикла next с оставшимся списком страниц: (next (cdr P)). Если числа есть, тогда берётся очередное простое число p. Если оно достаточно больше, or останавливается на истинном условии (fx> p l), и эта истина возвращается из обоих циклов. В противном случае, если p делит n, циклы возвращают ложь. А если p не делит n, простота n будет протестирована на следующей итерации loop очередным делителем на текущей странице: (loop (fx1+ i)).
Остаток кода: подсчёт количества чисел, записанных в список страниц, и основной цикл программы - должны быть понятны. Основной цикл реализует ту же логику, что и в prime-1.scm.
В первый вариант prime-2.scm прокралась ошибка (для тестирования возможностей компиляторов). Вместо того, чтобы вычислить границу для перебора простых делителей как n/2, то есть (fxquotient n 2), в оптимизаторском порыве я таки посчитал её, как (isqrt n), то есть, как ceil(sqrt(n)). От этого, естественно, объём вычислений и, соответственно, время работы резко сократились.
time scheme --script primes-3.scm 250000
22044
real 0m0.100s«Ну это уже ни в какие ворота!» - подумал я, увидев это, и решил, что с realloc или GCC что-то не так. Пытаясь исключить из тормозов realloc я переписал логику primes-2.scm на C. В процессе переписывания ошибку я заметил, но не пропадать же добру. primes-21.c - калька c primes-2.scm, с одной особенностью: страницы связываются в список не через отдельные cons-пары, а через указатели next, хранящиеся на самих страницах.
primes-21.c
Hidden text
// primes-21.c
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#define PAGE_CAPACITY 510
typedef struct page {
struct page *restrict next;
long cursor;
long numbers[PAGE_CAPACITY];
} Page;
Page *make_page(const long n) {
Page *const P = malloc(sizeof(Page));
P->next = NULL;
P->cursor = 1;
P->numbers[0] = n;
return P;
}
#define likely(X) __builtin_expect((X), 1) // __
Page *push(Page *restrict const P, const long n) {
if (likely((P->cursor < PAGE_CAPACITY))) {
P->numbers[P->cursor++] = n;
return P;
}
return P->next = make_page(n);
}
int is_prime(const Page *restrict const pages, const long n) {
const long l = n / 2;
for (const Page *restrict P = pages; P; P = P->next) {
for (int i = 0; i < P->cursor; i++) {
const long p = P->numbers[i];
if (p > l) return 1;
if (n % p == 0) return 0;
}
}
return 1;
}
long counting_free(const Page *restrict const pages) {
const Page *restrict P = pages;
long cnt = 0;
while (P) {
const Page *const next = P->next;
cnt += P->cursor;
free((void *)P);
P = next;
}
return cnt;
}
int main(int argc, char *argv[]) {
if (argc < 2) {
fprintf(stderr, "Specify limit\n");
return 1;
}
const long N = atol(argv[1]);
Page *restrict const primes = make_page(2);
Page *restrict P = primes;
for (long n = 3; n <= N; n += 2) {
if (is_prime(primes, n)) P = push(P, n);
}
printf("%ld\n", counting_free(primes));
return 0;
}
Итак, волнительный момент запуска.
$ gcc -O3 primes-2.c && time ./a.out 250000
22044
real 0m0.469s
$ gcc -O3 primes-21.c && time ./a.out 250000
22044
real 0m0.464s
$ time scheme --optimize-level 3 --script primes-2.scm 250000
22044
real 0m0.517s
«Во дела...» Код на C оказывается всего на 11% быстрее кода на Scheme. И оказывается, что работа realloc через механизмы виртуальной памяти не так уж сильно и помогает. Может быть, нужно считать больше, чтобы гарантированно выпасть из кэшей?
$ gcc -O3 primes-2.c && time ./a.out 2500000
183072
real 0m31.376s
$ gcc -O3 primes-21.c && time ./a.out 2500000
183072
real 0m31.323s
$ time scheme --optimize-level 3 --script primes-2.scm 2500000
183072
real 0m31.930s
2$ time pypy ../1/primes-1.py 2500000
183072
real 0m32.191sНемая сцена, в которой Жан-Люк Пикар молча прикладывает ладонь к лицу.
Итерация 3: Маловато будет.
Конечно, с границей для перебора делителей в n/2 много простых чисел не насчитаешь. Нужно перебирать до isqrt(n). Код программ итерации 3 только этим и отличается от кода программ итерации 2. Не буду приводить его здесь, он сохранён в gist.
С такими программами можно посчитать и до 25E+6.
$ g++ -O3 primes-3.cpp && time ./a.out 25000000
1565927
real 0m3.726s
$ gcc -O3 primes-3.c && time ./a.out 25000000
1565927
real 0m3.717s
$ gcc -O3 primes-31.c && time ./a.out 25000000
1565927
real 0m3.776s
$ time pypy primes-3.py 25000000
1565927
real 0m5.423s
$ time scheme --optimize-level 3 --script primes-3.scm 25000000
1565927
real 0m6.003s«Фьюх!..» Картина мира восстановлена. Программа на C в 1.45 раза быстрее программы на Python и в 1.6 раза быстрее программы на Scheme. Хотя, маловато, конечно, будет:
fixnum-ы в Scheme - не то же самое, что 64 битовые целые, арифметическая работа с fixnum-ами более замороченная;
структуа данных в программе на Scheme несколько сложнее;
Не полностью удовлетворённый результатом, я решил испытать Clang.
clang --version
Debian clang version 11.0.1-2
$ clang -O3 primes-3.c && time ./a.out 25000000
1565927
real 0m2.385s
$ clang -O3 primes-31.c && time ./a.out 25000000
1565927
real 0m2.388s
$ clang++ -O3 primes-3.cpp && time ./a.out 25000000
1565927
real 0m2.448s
$ clang primes-31.c && time ./a.out 25000000
1565927
real 0m4.667s
$ clang++ primes-3.cpp && time ./a.out 25000000
1565927
real 0m8.583s
$ clang primes-3.c && time ./a.out 25000000
1565927
real 0m4.590sКогда я увидел это, в глаза мне попала соринка то ли радости, то ли разочарования. Здесь требуется дальнейший анализ происходящего. Почему GCC выдаёт в 1.6 более медленный код, чем Clang? Почему неоптимизированная версия кода на C++ настолько хуже оптимизированной и при использовании Clang, и при использовании GCC? Но это требует вычитки получающегося машинного кода, на что у меня уже не оставалось времени.
Заключение
Scheme и Python языки динамические (PyPy работает только с программами, в которых он может статически вывести типы, но это программы на Python). Scheme совсем-совсем динамический, в этом языке даже нет обычных целых чисел фиксированной длинны: 32 битовых или 64 битовых. Fixnum-ы - значения в специальной кодировке, содержащей тип. Добавьте к этому сборку мусора, динамические контексты, замыкания, в которые превращается lambda, реализацию циклов (do - тоже макрос) через рекурсивные вызовы процедур (пусть и хвостовые).
С учётом этого, то, что PyPy и Chez Scheme умудряются выдавать код, который в данной конкретной числодробильной по своей природе задаче, работает всего лишь в 2.5 раза медленней, чем код, оптимизированный на 3-ем уровне системой LLVM (в которую вложено на пару порядков больше труда), удивляет.
Для себя из всего этого приключения я сделал вывод: имеет смысл писать прототипы алгоритмов на Python или Scheme, не только по той причине, что это проще и в написании, и в отладке, но и по той причине, что они позволяют получить некоторую базовую оценку производительности, к которой нужно стремиться. А то, ведь, есть риск, программируя на C++, получить, вроде как, естественную для задачи реализацию, которая будет раза в четыре медленней оптимальной.
Благодарю за внимание.
P.S. Если Вы хотите, чтобы я позапускал Ваши реализации primes на той самой виртуальной машине, чтобы можно было сравнить результаты, присылайте, пожалуйста, pull-реквесты сюда: https://github.com/mbakhterev/primes
