Показанным ниже кодом вы можете проверить на високосность год в интервале 0 ≤ y ≤ 102499 всего примерно тремя командами CPU:

bool is_leap_year_fast(uint32_t y) {
    return ((y * 1073750999) & 3221352463) <= 126976;
}

Как это работает? Ответ на удивление сложен. В статье я объясню процесс; в основном он связан с забавным битовым жонглированием. В конце мы обсудим применение этого кода на практике.

Вот, как обычно реализуется проверка на високосность:

bool is_leap_year(uint32_t y) {
    if ((y % 4) != 0) return false;
    if ((y % 100) != 0) return true;
    if ((y % 400) == 0) return true;
    return false;
}

Мы пользуемся пролептическим григорианским календарём, расширяющим григорианский календарь с момента начала его использования в 1582 году до нулевого года. Благодаря этому нам не нужно как-то по-другому обрабатывать даты до 1582 года. Для простоты мы не будем учитывать отрицательные года и используем беззнаковые величины.

Оптимизируем стандартное решение

Для начала реализуем простые трюки для ускорения, чтобы получить точку отсчёта. Я не знаю, кто конкретно их придумал — вероятно, эти трюки многократно изобретались заново.

Можно заменить (y % 100) != 0 на (y % 25) != 0: мы уже знаем, что y кратен 2², поэтому если он кратен 5², то кратен и 2² ⋅ 5² = 100. Аналогично, мы можем заменить (y % 400) == 0 на (y % 16) == 0: мы уже знаем, что y кратен 5², поэтому если он также кратен 2⁴, то кратен и 5² ⋅ 2⁴ = 400.

bool is_leap_year1(uint32_t y) {
    if ((y % 4) != 0) return false;
    if ((y % 25) != 0) return true;
    if ((y % 16) == 0) return true;
    return false;
}

Это полезно, потому что теперь мы можем заменить деление с остатком на 4 и 16 побитовым маскированием. Есть и ещё один трюк, хорошо известный разработчикам компиляторов, позволяющий избавиться от деления с остатком на 25. Скомпилировав (x % 25) != 0 при помощи gcc и выполнив трансляцию обратно на C, мы получим x * 3264175145 > 171798691. Так как умножение обычно имеет задержку 3 такта, а деление с остатком — не менее 20 тактов, это солидное улучшение. Я объясню в общих чертах, как это работает; подробности можно найти в

• В статье Faster remainder by direct computation: Applications to compilers and software libraries Дэниела Лемайра, Оуэна Кейсера и Натана Курца про уменьшение на константу делителя при делении с остатком;

• В Euclidean affine functions and their application to calendar algorithms Кассио Нери и Лоренца Шнайдера про календарные вычисления

• В Identifying leap years Дэвида Тёрнера о проверке високосности (в том числе с формальными доказательствами!).

Откуда берутся эти магические числа? У нас есть

2³² ⋅ 19/25 = 3264175144,96 (ровно).

То есть умножив на 3264175145, мы приблизительно вычисляем дробную часть умножения на (19/25). Если это было точное умножение, мы получим для кратных 25 чисел дробную часть, равную нулю. Однако если мы умножим на число, которое больше на 0,04, то погрешность может быть до 0,04 ⋅ (2³² - 1) = 171798691,8; отсюда и берётся второе магическое число.

Этот трюк не так хорошо работает для x % 100, где нам нужна ещё одна fixup-команда, так что уменьшение с y % 100 до y % 25 было очень полезным.

Теперь наша проверка на високосность выглядит так:

bool is_leap_year2(uint32_t y) {
    if ((y & 3) != 0) return false;
    if (y * 3264175145u > 171798691u) return true;
    if ((y & 15) == 0) return true;
    return false;
}

Стоит отметить, что современные компиляторы наподобие gcc и clang создают из is_leap_year1 что-то наподобие is_leap_year2, так что делать это в исходниках на C нет особого смысла, но в других языках это может оказаться полезным.

Обычно такой код компилируется в ассемблерный код с ветвлением. На практике вводимые данные обычно достаточно предсказуемы, поэтому это необязательно плохо. Если мы хотим избежать снижения производительности из-за ошибок прогнозирования ветвления, то ценой замедления наилучшего случая можем немного изменить структуру и получить код без ветвления (и одновременно хорошего кандидата для соревнований в гольф-кодинге):

bool is_leap_year3(uint32_t y) {
    return !(y & ((y % 25) ? 3 : 15));
}

Если вы хотите узнать о других способах ускорения календарных вычислений, то изучите Optimizing with Novel Calendrical Algorithms Джейкоба Пратта.

Находим решение с битовым жонглированием

Можно ли усовершенствовать вычисление високосности, отказавшись от корректности всех входных данных? В конце концов, обычно нас не волнует, будет ли високосным 3584536493 год; и в самом деле, Python, C# и Go поддерживают года с 0 (или 1) до 9999 (в котором смещение относительно времён года уже будет больше четырёх дней). Я подумал, что если существует более короткое решение, то оно бы выглядело как какое-то странное хэширование с магическими константами, поэтому решил попробовать малые формы и подобрать константы перебором. Вид (y * f) <= t кажется полезным, но пока недостаточно мощным. Одним из моих кандидатов стало добавление маски: ((y * f) & m) <= t. Теперь у нас есть для угадывания 96 бита, одним лишь брутфорсом это не решить. Давайте воспользуемся z3 — солвером, поддерживающим ограничения битовых векторов, что идеально нам подходит:

import z3

BITS = 32
f, m, t, y = z3.BitVecs('f m t y', BITS)

def target(y):
    return z3.And((y & 3) == 0, z3.Or(z3.URem(y, 25) != 0, (y & 15) == 0))

def candidate(x):
    return z3.ULE((x * f) & m, t)

solver = z3.Solver()
solver.add(z3.ForAll(y, z3.Implies(z3.ULE(y, 400),
                                   candidate(y) == target(y))))

if solver.check() == z3.sat:
    print(f'found solution: {solver.model()}')
else:
    print('no solution found')

За несколько секунд он обнаружил константы, позволяющие получать корректный результат для нетривиального интервала значений годов. Расширяя этот интервал, я в конечном итоге примерно за полчаса вычислений обнаружил константы, дающие корректный результат с года 0 по год 102499, и доказал, что для 32 битов это оптимум:

bool is_leap_year_fast(uint32_t y) {
    const uint32_t f = 1073750999u;
    const uint32_t m = 3221352463u;
    const uint32_t t = 126976u;
    return ((y * f) & m) <= t;
}

Объяснение

Как это работает? Кажется неожиданным, почти волшебным, что мы можем уместить все эти вычисления всего в три команды. Однако описанное выше даёт нам основные инструменты для понимания этого.

Вот наши константы в двоичном виде; четыре релевантных интервала битов обозначены буквами:

Давайте сначала разберёмся, что делает маскирование при помощи m и сравнение с t для произведения p := y ⋅ f. В блоке A биты t равны 0, а значит, если какой-то из битов в A ненулевой в p, результатом будет false. В противном случае релевантным становится блок B. Здесь все биты в t равны 1, поэтому результат равен true, если какой-то из битов B имеет в p ненулевое значение. В противном случае для блока C мы требуем, чтобы все биты в p были нулевыми. Благодаря этому куча сравнений интервалов битов объединена в единственный <=.

То есть мы можем переписать is_leap_year_fast следующим образом:

bool is_leap_year_fast2(uint32_t y) {
    uint32_t p = y * 1073750999u;
    const uint32_t A = 0b11000000000000000000000000000000;
    const uint32_t B = 0b00000000000000011111000000000000;
    const uint32_t C = 0b00000000000000000000000000001111;
    if ((p & A) != 0) return false;
    if ((p & B) != B) return true;
    if ((p & C) == 0) return true;
    return false;
}

Это подозрительно похоже на is_leap_year2! И в самом деле, три условия имеют точно такое же предназначение. Мы покажем, что

1. (p & A) != 0 срабатывает, когда (y % 4) != 0;

2. (p & B) != B срабатывает, когда (y % 100) != 0;

3. (p & C) == 0 срабатывает, когда (y % 16) == 0 (и, следовательно, (y % 400) == 0, потому что мы уже знаем, что y кратно 25).

Два простых случая: (1) и (3)

(1): Бит 1 блока A в f воссоздаёт два младших бита y в p в блоке A. Их нельзя спутать с результатом умножения на биты в D: максимальное значение, которое мы можем получить — это 102499 ⋅ (f & D) = 940428325, состоящее всего из 30 битов. Таким образом, проверка A на нулевое значение в p эквивалентно проверке, равен ли нулю остаток от деления y на 4.

(3): Проверка того, что ни один из младших 4 битов не задан в p — это проверка того, равен ли нулю остаток от деления p на 16. Однако на самом деле мы хотим проверить y. Это не проблема: достаточно взглянуть только на младшие 4 бита f, а f здесь 1111₂ = 15. Умножение на 15 = 3 ⋅ 5 не добавляет нового делителя 2, поэтому состояние делимости на 16 не меняется.

Любопытный случай: (2)

Теперь давайте попробуем выяснить, для каких чисел p & B ≠ B. Для этого бит 1 в f & A не играет роли, поэтому рассмотрим биты в f & D. Они равны 10001111010111₂ = 9175. Давайте проверим, какие числа пройдут тест.

>>> B = 0b00000000000000011111000000000000
>>> s = [y for y in range(5000) if ((y * 9175) & B) == B]
>>> for i in range(0, len(s), 16): print(*(f'{n:4d}' for n in s[i:i+16]))
  14   57   71  100  114  157  171  200  214  257  271  300  314  357  371  400
 414  457  471  500  514  557  571  600  614  657  671  700  714  757  771  800
 814  857  871  900  914  957  971 1000 1014 1057 1071 1100 1114 1157 1171 1200
1214 1257 1271 1300 1314 1357 1371 1400 1414 1457 1471 1500 1514 1557 1571 1600
1614 1657 1671 1700 1714 1757 1771 1800 1814 1857 1871 1900 1914 1957 1971 2000
2014 2057 2071 2100 2114 2157 2171 2200 2214 2257 2271 2300 2314 2357 2371 2400
2414 2457 2471 2500 2514 2557 2571 2600 2614 2657 2671 2700 2714 2757 2771 2800
2814 2857 2871 2900 2914 2957 2971 3000 3014 3057 3071 3100 3114 3157 3171 3200
3214 3257 3271 3300 3314 3357 3371 3400 3414 3457 3471 3500 3514 3557 3571 3600
3614 3657 3671 3700 3714 3757 3771 3800 3814 3857 3871 3900 3914 3957 3971 4000
4014 4057 4071 4100 4114 4157 4200 4214 4257 4300 4314 4357 4400 4414 4457 4500
4514 4557 4600 4614 4657 4700 4714 4757 4800 4814 4857 4900 4914 4957

Здесь есть числа, кратные 100, как мы и хотели, но и куча других чисел. Это не проблема, если ни одно из них не окажется кратным 4, потому что мы уже отфильтровали их на предыдущем этапе. Кроме того, отсутствует 0, но это тоже не проблема, потому что 0 также кратен 400.

Давайте попробуем разобраться в паттерне. На первый взгляд, он выглядит очень простым: у нас есть *14, *57, *71, и *00. Однако начиная с 4171 числа *71 пропадают (вы заметили?). Позже появляются новые паттерны. Давайте проанализируем это. Результат работы программы на Python

def test(y):
    B = 126976
    return ((y * 9175) & B) == B

active = set()
for y in range(120000):
    r = y % 100
    if test(y):
        if r not in active:
            print(f'{y:6}: started *{r:02}')
            active.add(r)
    else:
        if r in active:
            print(f'{y:6}: stopped *{r:02}')
            active.remove(r)

будет таким:

    14: started *14
    57: started *57
    71: started *71
   100: started *00
  4171: stopped *71
 32843: started *43
 36914: stopped *14
 65586: started *86
 69657: stopped *57
 98329: started *29
102500: stopped *00

То есть начиная с 102500 мы больше не отлавливаем числа, кратные 100, и именно поэтому 102499 — последнее число, для которого is_leap_year_fast возвращает корректный результат. Также мы видим, что ниже него нет ни одного числа, кратного 4, за исключением кратных 100 (удобно, что мы можем проверить это, зная только два последних десятичных разряда). Если мы доверимся этой проверке перебором, то доказательство условия (2) на этом завершается; но давайте продолжим, чтобы лучше понять, почему получаются именно такие числа.

Для начала давайте разберёмся, почему мы получаем значения, кратные 100. Делитель 9175 близок к кратному 1/100 в 17-битном представлении с фиксированной запятой:

2¹⁷ ⋅ 7/100 = 9175,04 (ровно).

Умножив число, кратное 100, на 9175,04, мы получим целое число (кратное 7) в битах 17 и выше, а биты ниже 17 будут нулевыми. Пример:

9175,04 ⋅ 500 = 10001100000000000000000₂, где 100011₂ = 35 = 5 ⋅ 7.

Умножив число, кратное 100, на 9175, мы получим немного меньше:

9175 ⋅ 500 = 10001100000000000000000₂ − 500 ⋅ 0,04 = 10001011111111111101100₂.

В общем случае, вычитание небольшого значения из числа, заканчивающегося на кучу нулей, создаёт число, заканчивающееся на кучу единиц, за исключением самого конца. Здесь мы проверяем 5 битов в B. Для y, кратного 100, они все гарантированно будут единицами, если только накопленная ошибка не достигнет конца B, что произойдёт только после y = 2¹² / 0,04 = 102400, что нам подходит.

Откуда же берутся другие числа наподобие 14, 57 и 71? Давайте взглянем на это под другим углом. У нас есть 9175 = 2¹⁷ ⋅ 0,06999969482421875 (ровно) и B = 2¹⁷ ⋅ 0,96875, поэтому

p & B= B

⇔{y ⋅ 0,06999969482421875}≥ 0,96875, где {x} — дробная часть x

⇔6,999969482421875y mod 100≥ 96,875

Это ещё один способ понять, почему принимаются числа, кратные 100: для них 7y mod 100 равно 0, поэтому 6,999969482421875y mod 100 оказывается чуть меньше 100, и падает ниже 96,875 только после y = (100 − 96,875) / (7 − 6,999969482421875) = 102400.

Чтобы понять другие числа, встречающиеся в нашей последовательности, давайте сначала рассмотрим, какими будут решения, если бы в этом неравенстве было ровно 7:

7y mod 100≥ 96,875

⇔ 7y mod 100∈ {97, 98, 99}.

Чтобы найти решения этого, нам сначала потребуется число, обратное по модулю 7 modulo 100, то есть число x такое, что 7x mod 100 = 1. Мы можем вычислить его при помощи расширенного алгоритма Евклида или просто при помощи какого-нибудь онлайн-калькулятора, который сообщит нам, что результат равен 43. Тогда решениями будут 43 ⋅ 97 (mod 100), 43 ⋅ 98 (mod 100) и 43 ⋅ 99 (mod 100), то есть, соответственно, 71, 14 и 57 (mod 100). Это объясняет, почему мы сначала встречаем числа вида *14, *57 и *71. Это также объясняет, почему мы перестаём встречать, например, *71 после 4071: хотя 7 ⋅ 4171 = 29197, мы получаем 6,999969482421875 ⋅ 4171 = 29196,872711181640625, что (modulo 100) меньше, чем 96,875. Аналогично, мы встречаем 32843, потому что накопившаяся погрешность (7 − 6,999969482421875) ⋅ 32843 = 1,002288818359375 превышает единицу. Приложив ещё немного усилий, мы можем вручную воссоздать результат приведённой выше программы на Python и убедиться, что ни одно из этих чисел не кратно 4.

Расширение до других битовых ширин

Теперь, когда мы знаем, как работает этот трюк, можно попробовать подобрать параметры для других битовых ширин. Изменятся местоположение блока B и числитель знаменателя 100 в f & D.

uint64_t score(uint64_t f, uint64_t m, uint64_t t) {
      for (uint64_t y = 0; ; y++)
          if ((((y * f) & m) <= t) != is_leap_year(y))
              return y;
  }
  
  int main() {
      uint64_t best_score = 0;
      for (int k = 0; k < BITS; k++) {
          for (int k2 = 0; k2 < k; k2++) {
              uint64_t t = (1ULL << k) - (1ULL << k2);
              uint64_t m = (0b11ULL << (BITS - 2)) | t | 0b1111;
              for (int n = 0; n < 100; n++) {
                  uint64_t f = (0b01ULL << (BITS - 2)) | (((1ULL << k) * n) / 100);
                  uint64_t new_score = score(f, m, t);
                  if (new_score > best_score) {
                      printf("%llu %llu %llu: %llu (%d %d %d)\n",
                             f, m, t, new_score, k, k - k2, n);
                      best_score = new_score;
                  }
              }
          }
      }
      return 0;
  }

При BITS = 64 мы примерно за 7 минут находим f = 4611686019114582671, m = 13835058121854156815, t = 66571993088, которые корректны вплоть до y = 5965232499. Это здорово, потому что 5965232499 > 2³², а значит, таким вариантом кода можно протестировать любой 32-битный год.

Какого максимального года мы можем достичь с 64 битами? Возможно, есть другие константы, которые работают ещё лучше? Я не могу сходу найти способ доказать это, поэтому попросил заняться этим других людей, создав пост о задаче в Code Golf StackExchange. И спустя всего один час пользователь ovs опубликовал очень хороший результат, а два дня спустя пользователь Exalted Toast выложил доказательство того, что 5965232499 и в самом деле — наилучший возможный интервал для 64 битов, тоже воспользовавшись солвером z3.

Бенчмарк

Получить чёткие бенчмарки в этом случае сложно, потому что функция выполняется очень малое количество; более того, время исполнения версий с ветвлением зависит от паттернов ввода. Попробуем два крайних случая: всегда 2025 год и полностью случайные годы. Ниже приведены результаты бенчмарка на i7-8700K (Coffee Lake, 4.7 GHz), скомпилированного с g++ -O3 -fno-tree-vectorize:

Вот некоторые из странностей:

• При случайных значениях is_leap_year2 чуть медленнее is_leap_year. Это неожиданно, потому что для y % 100 требуется на одну команду больше, чем трюку в is_leap_year2.

• is_leap_year3 чуть медленнее для случайных данных, чем для фиксированного значения. Это неожиданно, потому что функция не выполняет никакого ветвления.

У меня нет никаких других объяснений этому, кроме как то, что создание бенчмарков — сложная задача.

Для случайных данных новая функция is_leap_year_fast в 3,8 раза быстрее, чем стандартная реализация, а для полностью прогнозируемого ввода она примерно на 6% медленнее. В целом, результаты кажутся вполне стабильными.

Заключение

Стоит ли оно того? Нужно ли нам заменять, например, реализацию datetime CPython на этот трюк? Ответ зависит от обстоятельств. На практике, чаще всего будут запрашивать проверку текущего года, или, по крайней мере, запросы будут достаточно предсказуемы; в таком случае особого преимущества мы не получим. Чтобы изменения оправдали себя, в идеале нам нужен бенчмарк с реалистичными данными, который использует в качестве подпрограммы проверку года на високосность, а не простой микробенчмарк. С удовольствием услышал бы о результатах подобных измерений!