Julian map と packed table による高速な civil_from_days
Unix epochからの通算日を年月日へ戻す処理は、一般に civil_from_days と呼ばれます。前の記事で扱った days_from_civil の逆変換です。
civil date -- days_from_civil --> epoch day
civil date <-- civil_from_days --- epoch dayこの記事では出力年を 1..9999 に限定し、Ben JoffeのJulian mapによる年復元と、366要素のpacked month/day tableを組み合わせます。Bunではさらに、二つの定数除算を上方丸めしたbinary64逆数との乗算へ置き換えます。
Apple Silicon上の今回の測定では、実行時間が通常のHoward Hinnant実装の約13%（Bun）、26%（Rust）、30%（C）になりました。削減率ではそれぞれ約87%、74%、70%です。
前提
入力範囲は、暦年 1..9999 に対応する次のUnix epoch dayです。
-719162 <= epochDay <= 2932896両端はそれぞれ 0001-01-01 と 9999-12-31 です。
TypeScript実装
Bun向けのホットパスは、年月日を一つの整数へpackして返します。
packed = year * 512 + month * 32 + day月は 1..12、日は 1..31 なので、下位9 bitへ月日を格納できます。
const INV_146097 = (1 / 146_097) * (1 + Number.EPSILON);
const INV_1461 = (1 / 1_461) * (1 + Number.EPSILON);
const MONTH_DAY = new Uint16Array(366);
for (let dayOfMarchYear = 0; dayOfMarchYear <= 365; dayOfMarchYear++) {
const n = 2_141 * dayOfMarchYear + 197_913;
const marchMonth = n >>> 16;
const day = Math.floor((n & 65_535) / 2_141) + 1;
const janOrFeb = dayOfMarchYear >= 306;
const month = janOrFeb ? marchMonth - 12 : marchMonth;
MONTH_DAY[dayOfMarchYear] =
(Number(janOrFeb) << 9) | (month << 5) | day;
}
/**
* Returns year * 512 + month * 32 + day.
* Preconditions: -719162 <= epochDay <= 2932896.
*/
export function civilFromDaysPacked(epochDay: number): number {
const q = ((epochDay << 2) + 2_877_875) | 0;
const century = (q * INV_146097) | 0;
const julian =
(q + ((century - (century >> 2)) << 2)) | 0;
const year = (julian * INV_1461) | 0;
const remainder =
(julian - Math.imul(year, 1_461)) | 0;
return (
(year << 9) + MONTH_DAY[remainder >>> 2]
) | 0;
}通常のオブジェクトが必要なら、ホットループの外で展開できます。
export function unpackCivilDate(packed: number) {
return {
year: Math.floor(packed / 512),
month: (packed >>> 5) & 15,
day: packed & 31,
};
}civilFromDaysPacked 自体にオブジェクト生成を含めると、カレンダー計算より割り当てコストを測る比率が大きくなるため、ベンチマークではpackした整数を返しています。
Rust実装
Rustではtableをコンパイル時に生成できます。定数除算はそのまま記述し、最適化をコンパイラに任せます。
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct CivilDate {
pub year: i32,
pub month: u32,
pub day: u32,
}
const fn build_month_day_table() -> [u16; 366] {
let mut table = [0u16; 366];
let mut r = 0usize;
while r <= 365 {
let n = 2_141 * r as u32 + 197_913;
let march_month = n >> 16;
let day = (n & 65_535) / 2_141 + 1;
let month = if r >= 306 {
march_month - 12
} else {
march_month
};
let year_bump = u16::from(r >= 306);
table[r] = (year_bump << 9) | ((month << 5) | day) as u16;
r += 1;
}
table
}
const MONTH_DAY: [u16; 366] = build_month_day_table();
#[inline]
pub fn civil_from_days(epoch_day: i32) -> CivilDate {
let q = 4 * (epoch_day + 719_468) as u32 + 3;
let century = q / 146_097;
let julian = q + century * 3 + (century & 3);
let year = julian / 1_461;
let day_of_march_year = julian % 1_461 / 4;
let packed = MONTH_DAY[day_of_march_year as usize] as u32;
CivilDate {
year: (year + (packed >> 9)) as i32,
month: (packed >> 5) & 15,
day: packed & 31,
}
}C実装
Cでは次の初期化をプログラム開始時に一度実行します。本番用途では、生成結果を static const uint16_t[366] として埋め込めば初期化も不要です。
#include <stdint.h>
typedef struct {
int32_t year;
uint32_t month;
uint32_t day;
} civil_date_t;
static uint16_t month_day[366];
static void init_month_day_table(void) {
for (uint32_t r = 0; r <= 365; ++r) {
const uint32_t n = 2141u * r + 197913u;
const uint32_t march_month = n >> 16;
const uint32_t day = (n & 65535u) / 2141u + 1u;
const uint32_t month = r >= 306u
? march_month - 12u
: march_month;
month_day[r] = (uint16_t)(
((r >= 306u) << 9) | (month << 5) | day
);
}
}
static inline civil_date_t civil_from_days(int32_t epoch_day) {
const uint32_t q = 4u * (uint32_t)(epoch_day + 719468) + 3u;
const uint32_t century = q / 146097u;
const uint32_t julian = q + century * 3u + (century & 3u);
const uint32_t year = julian / 1461u;
const uint32_t day_of_march_year = julian % 1461u / 4u;
const uint32_t packed = month_day[day_of_march_year];
return (civil_date_t) {
(int32_t)(year + (packed >> 9)),
(packed >> 5) & 15u,
packed & 31u
};
}通常の civil_from_days
Howard Hinnantの基準実装をTypeScriptで書くと次のようになります。
function civilFromDaysHinnant(epochDay: number) {
const z = epochDay + 719_468;
const era = Math.floor(z / 146_097);
const dayOfEra = z - era * 146_097;
const yearOfEra = Math.floor(
(dayOfEra -
Math.floor(dayOfEra / 1_460) +
Math.floor(dayOfEra / 36_524) -
Math.floor(dayOfEra / 146_096)) /
365,
);
const year = yearOfEra + era * 400;
const dayOfYear =
dayOfEra -
(365 * yearOfEra +
Math.floor(yearOfEra / 4) -
Math.floor(yearOfEra / 100));
const marchMonth = Math.floor((5 * dayOfYear + 2) / 153);
const day =
dayOfYear - Math.floor((153 * marchMonth + 2) / 5) + 1;
const month = marchMonth + (marchMonth < 10 ? 3 : -9);
return {
year: year + (month <= 2 ? 1 : 0),
month,
day,
};
}これは広い年範囲を扱える、明快で移植性の高い実装です。今回の限定版は、この汎用性を固定範囲での速度と交換します。
Julian mapで年を復元する
今回の年復元はBen JoffeのJulian mapを利用します。
q = 4 * (epochDay + 719468) + 3
century = floor(q / 146097)
julian = q + 3 * century + (century & 3)
year = floor(julian / 1461)
rem = floor((julian - 1461 * year) / 4)Gregorian calendarの400年周期を補正した後、一時的にJulian calendarの4年周期へ写すことで、年復元の依存鎖を短くします。
元のJoffe式 q - (century & ~3) + 4 * century は、century = (century & ~3) + (century & 3) を使うと、上の短い式へ整理できます。
1461 = 4 * 365 + 1 なので、julian / 1461 がMarch-basedの年を返し、余りを4で割った値がMarch 1から数えた年内日になります。
0 <= rem <= 3651月と2月はMarch-based yearの末尾にあります。rem >= 306 の年補正をtableのbit 9へ格納することで、実行時の比較と加算も除去します。
月日を366要素のtableから取得する
通常は年内日から月日を復元するために、乗算、シフト、除算、1月・2月の補正を行います。しかし rem は 0..365 の366通りしかありません。
各値に対応する年補正、月、日を次の10-bit整数へ格納します。
(yearBump << 9) | (month << 5) | day取得は次だけです。
const packed = MONTH_DAY[dayOfMarchYear];
const year = marchYear + (packed >>> 9);
const month = (packed >>> 5) & 15;
const day = packed & 31;tableの大きさは、
366 * 2 = 732 bytesです。一般的なPCのL1 cacheには十分小さい一方、cold cache、組み込み環境、コードサイズを重視する用途では算術版が有利な可能性があります。
Bunで除算を上方丸め逆数へ置き換える
Bunでは /146097 と /1461 より、あらかじめ求めた逆数との乗算が高速でした。
const INV_146097 = (1 / 146_097) * (1 + Number.EPSILON);
const INV_1461 = (1 / 1_461) * (1 + Number.EPSILON);実際のbinary64値は次のとおりです。
INV_146097 = 0.0000068447675174712704
bits = 0x3edcb5835e647c33
INV_1461 = 0.00068446269678302542
bits = 0x3f466db072f2284e通常の最も近いbinary64逆数では、割り切れる入力の積が整数よりわずかに小さくなり、Math.floor が1小さい値を返す場合があります。1 + Number.EPSILON を掛け、真の逆数よりわずかに大きい隣接値を使うことでこれを避けます。
どちらも、
INV_d * d - 1 = 2^-52です。対象範囲の最大商は約100と10000なので、上方誤差はそれぞれ約 2.3e-14 と 2.3e-12 以下です。一方、割り切れない整数入力が次の商境界まで持つ距離は少なくとも 1/146097 と 1/1461 です。したがって、上方誤差によって次の整数へ到達することはありません。
この議論は指定範囲とIEEE-754 binary64に依存します。範囲を広げる場合は再証明と再検証が必要です。
また、二つの積は非負で2³¹未満の整数へ切り捨てられるため、Bun版では Math.floor を32-bit整数化の | 0 に置き換えられます。
const century = (q * INV_146097) | 0;
const year = (julian * INV_1461) | 0;最終版では、q の定数部分も畳み込み、シフト、積、戻り値をint32ドメインへ閉じています。
const q = ((epochDay << 2) + 2_877_875) | 0;
const remainder =
(julian - Math.imul(year, 1_461)) | 0;また、非負の century について、
3c + (c & 3) = 4 * (c - floor(c / 4))なので、Bun版では次の別形を使います。
const julian =
(q + ((century - (century >> 2)) << 2)) | 0;これは今回のBunでは有効でしたが、最適化されたC/Rustではコンパイラが元の式をすでに同等の命令列へ変換できます。一般にこちらが速いと主張するものではありません。
追加tableなしの4-wide SIMD版
複数の日付をまとめて変換できる場合、732-byte tableをさらに大きくする代わりに、月日復元を算術式へ戻して4入力をNEONで並列処理できます。
January/February補正は、packed表現では次の一回の加算になります。
year += 1, month -= 12
packed += 512 - 12 * 32
= 128Clangの4要素vector型で中心部分を書くと次の形です。
typedef uint32_t u32x4 __attribute__((ext_vector_type(4)));
static inline u32x4 civil_from_days_4(u32x4 epoch_day) {
const u32x4 q = epoch_day * 4u + 2877875u;
const u32x4 century = q / 146097u;
const u32x4 julian = q + century * 3u + (century & 3u);
const u32x4 year = julian / 1461u;
const u32x4 rem = (julian - year * 1461u) >> 2;
const u32x4 n = rem * 2141u + 197913u;
const u32x4 march_month = n >> 16;
const u32x4 day = (n & 65535u) / 2141u + 1u;
const u32x4 jan_feb = (u32x4)(rem >= 306u);
return year * 512u + march_month * 32u + day
+ (jan_feb & 128u);
}Apple Silicon／Clang -O3 -march=native ではNEON命令へ展開されました。同じランダム入力配列を処理し、scalar側の自動vectorizeを無効にしたbatch比較は次のとおりです。
| 算術版 | ns/件 | scalar比 |
|---|---|---|
| scalar | 2.694 | 1.00倍 |
| 4-wide SIMD | 0.713 | 3.78倍 |
全3,652,059 epoch dayについてHinnant実装と一致しました。これはARM64向けbatch APIのスループットであり、1件を返すscalar APIやpacked table版のレイテンシと同一の比較ではありません。
32-bit整数の境界
入力範囲では、
306 <= epochDay + 719468 <= 3652364
1227 <= q <= 14609459
0 <= century <= 99
0 <= julian < 14610000
0 <= year <= 9999
0 <= dayOfMarchYear <= 365となり、整数部分は符号付き32-bitへ十分収まります。year bumpを含む MONTH_DAY の各値は10 bit、table全体は732 bytesです。
全数検証
指定範囲のすべてのepoch dayについて、Hinnant実装と比較しました。
2932896 - (-719162) + 1
= 3,652,059 inputs検証対象はHinnant、Neri–Schneider、Ben Joffe、Bunの上方逆数版、packed table版です。TypeScript、Rust、Cの最終実装は全件一致しました。
ベンチマーク
Apple Silicon ARM64上で、決定的な疑似乱数で生成した 2^20 件のepoch dayを24回処理しました。入力位置はラウンドごとに回転させています。5回測定した最小値です。
Bun 1.4.0
rustc 1.96.0 -C opt-level=3 -C target-cpu=native
Apple clang 21.0.0 -O3 -march=native今回再測定した結果は次のとおりでした。
| 環境 | 通常のHinnant | 今回の限定版 | 実行時間削減率 | スループット |
|---|---|---|---|---|
| Bun | 37.088 ns/件 | 4.830 ns/件 | 87.0% | 7.68倍 |
| Rust | 5.734 ns/件 | 1.500 ns/件 | 73.8% | 3.82倍 |
| C | 5.315 ns/件 | 1.603 ns/件 | 69.8% | 3.32倍 |
Bunでは、year bumpをtableへ移しつつdouble演算を残した版が6.091 ns/件、int32ドメインへ閉じた最終版が4.830 ns/件でした。year bumpのtable化に加え、定数畳み込み、Math.imul、シフトを使ってdoubleとint32の往復を減らした効果が出ています。
RustとCでは、ソース上の定数除算をコンパイラが乗算とシフトへ変換できます。そのため、明示的な逆数より通常の整数除算表記が適切でした。packed tableによる月日復元の依存鎖削減が大きく効いています。
これらは特定のCPU、コンパイラ、JIT、入力分布での結果です。別世代のApple Silicon、x86-64、Node.js/V8、ブラウザ、別バージョンのRustやClangでは必ず再測定してください。
先行研究と今回の位置づけ
この実装の年復元はBen JoffeのJulian mapを利用しています。Gregorian calendar変換を乗算とシフトへ変換する方法は、NeriとSchneiderがEuclidean affine functionとして体系化しています。Howard Hinnantの実装は基準となる汎用アルゴリズムです。
したがって、Julian map、年復元、定数除算のstrength reduction、年内日tableのいずれも、単独の新規技法として主張するものではありません。
今回の構成は次の限定特殊化です。
- 出力年を
1..9999に限定する - Ben JoffeのJulian mapで年を復元する
0..365のMarch-based年内日を、year bumpも含む732-byteのpacked tableへ写す- Bunでは範囲証明した上方丸めbinary64逆数を使う
- CとRustでは定数除算のstrength reductionをコンパイラに任せる
調査した公開実装の範囲では、この固定範囲、Bunの上方丸め逆数、Joffe式、packed month/day tableを組み合わせた同一構成は確認できませんでした。ただし、あらゆるCPU、言語、入力範囲で最速または最適であることを意味しません。
正確な表現は次のとおりです。
Apple Silicon上の今回のBun、Rust、Cベンチマークにおいて、比較したHinnant、Neri–Schneider、Ben Joffe実装より高速だった。
まとめ
年 1..9999 に対応するepoch dayへ範囲を限定すると、civil_from_days は次の構成で高速化できました。
- Julian mapでGregorianの400年補正を行う
- 1461日のJulian 4年周期から年とMarch-based年内日を復元する
- 366要素のpacked tableからyear bump、月、日を一度に取得する
- Bunでは二つの定数除算を上方丸め逆数との乗算へ置き換える
今回のApple Silicon測定では、実行時間が通常のHinnant実装の約13%（Bun）、26%（Rust）、30%（C）になりました。さらに、追加tableを使わないCの4-wide SIMD算術版は、scalar算術版の約3.78倍のbatch throughputでした。
汎用コードではHinnant実装、tableを避けるならNeri–SchneiderまたはJoffeの算術版が適しています。入力範囲が保証され、732-byte tableが許容でき、変換がホットパスにある場合に今回の限定版が候補になります。
参考資料
- Howard Hinnant, chrono-Compatible Low-Level Date Algorithms
- Cassio Neri and Lorenz Schneider, Euclidean Affine Functions and Applications to Calendar Algorithms
- Cassio Neri, EAF supplementary implementation
- Ben Joffe, The Julian Map
- Ben Joffe, A Very Fast 64-Bit Date Algorithm
- Ben Joffe, fast-date-benchmarks
- マジックナンバー化による高速な
days_from_civil