マジックナンバー化による高速な曜日計算
曜日を求めるアルゴリズムとして、Tomohiko Sakamotoの手法がよく知られています。
この記事では、対象年を1年から9999年に限定し、Sakamotoの計算に含まれる定数除算とmod 7を、乗算・シフト・ビット抽出に置き換えます。
Apple Silicon上でCのインライン・バッチ処理を測定したところ、通常のSakamoto実装に対して約30%の高速化を確認しました。
最初に完成した実装を紹介し、その後で定数5243と74899が機能する理由を説明します。
実装例
TypeScript
/** Returns 0=Sunday, ..., 6=Saturday.
*
* Preconditions:
* - year: 1..9999
* - month: 1..12
* - day: 1..31
*/
const SCALED_MONTH_TERM = new Uint32Array([
0, 224_697, 149_798, 374_495, 0, 224_697, 374_495, 74_899, 299_596, 449_394, 149_798, 299_596,
]);
export function weekdayBoundedDistributed(year: number, month: number, day: number): number {
year -= Number(month < 3);
const product = Math.imul(year, 5_243);
const base = year + (year >>> 2) - (product >>> 19) + (product >>> 21) + day;
return ((Math.imul(base, 74_899) + SCALED_MONTH_TERM[month - 1]) >>> 16) & 7;
}TypeScriptでは、32-bit整数乗算を明示するためにMath.imulを使用しています。また、右シフトには符号なし右シフトの>>>を使用しています。
Rust
/// Returns 0=Sunday, ..., 6=Saturday.
///
/// Preconditions:
/// - year: 1..9999
/// - month: 1..12
/// - day: 1..31
#[inline]
pub fn weekday_bounded_distributed(
mut year: u32,
month: u32,
day: u32,
) -> u32 {
const SCALED_MONTH_TERM: [u32; 12] = [
0, 224_697, 149_798, 374_495,
0, 224_697, 374_495, 74_899,
299_596, 449_394, 149_798, 299_596,
];
year -= (month < 3) as u32;
let product = year * 5_243;
let base = year
+ (year >> 2)
- (product >> 19)
+ (product >> 21)
+ day;
(
(
base * 74_899 +
SCALED_MONTH_TERM[(month - 1) as usize]
) >> 16
) & 7
}C
#include <stdint.h>
/*
* Returns 0=Sunday, ..., 6=Saturday.
*
* Preconditions:
* - year: 1..9999
* - month: 1..12
* - day: 1..31
*/
static inline uint32_t weekday_bounded_distributed(
uint32_t year,
uint32_t month,
uint32_t day
) {
static const uint32_t scaled_month_term[12] = {
0u, 224697u, 149798u, 374495u,
0u, 224697u, 374495u, 74899u,
299596u, 449394u, 149798u, 299596u
};
year -= month < 3;
const uint32_t product = year * 5243u;
const uint32_t base = year
+ (year >> 2)
- (product >> 19)
+ (product >> 21)
+ day;
return (
(
base * 74899u +
scaled_month_term[month - 1]
) >> 16
) & 7u;
}通常のSakamotoの手法
Sakamotoの手法では、次の月補正値を使用します。
0, 3, 2, 5, 0, 3, 5, 1, 4, 6, 2, 4Cで書くと、基本形は次のようになります。
static uint32_t weekday_sakamoto(
uint32_t year,
uint32_t month,
uint32_t day
) {
static const uint8_t month_term[12] = {
0, 3, 2, 5, 0, 3,
5, 1, 4, 6, 2, 4
};
year -= month < 3;
return (
year
+ year / 4
- year / 100
+ year / 400
+ month_term[month - 1]
+ day
) % 7;
}1月と2月を前年の13月、14月に相当するものとして扱うため、month < 3の場合に年を1減らします。
その後、
year
+ floor(year / 4)
- floor(year / 100)
+ floor(year / 400)
+ month term
+ dayを計算し、7で割った余りを曜日とします。
Sakamotoの手法については、次の記事でも実装と考え方が紹介されています。
マジックナンバー化
今回の検討は、Falk Hüffner氏による次の記事をきっかけに始めました。
この記事では、入力範囲を限定することで、うるう年判定を乗算、マスク、比較による短い式へ変換しています。
同じように、曜日計算でも入力範囲を限定すれば、Sakamoto式の定数除算や剰余演算を短い整数演算へ変換できるのではないかと考えました。
今回の実装では、年を1以上9999以下に限定しています。
1月と2月について年を1減らした後でも、年の値は次の範囲に収まります。
0 <= year <= 9999この限定範囲を利用して、次の二つの変換を行います。
const uint32_t product = year * 5243u;
product >> 19; // year / 100
product >> 21; // year / 400さらに、Sakamoto式の中間値をvalueとすると、次の式でvalue % 7を求められます。
((value * 74899u) >> 16) & 7u今回のポイントは、次の二つの定数です。
5243
74899ただし、5243を使った0以上9999以下の整数の100除算は既知の最適化(strength reduction)です。
今回の実装では、同じ積year * 5243からyear / 100とyear / 400の両方を取り出しています。こちらについては後ほど説明します。
月補正値の分配
最初に作成した実装では、通常の月補正値を加えてから、全体を74899倍していました。
value =
year_term
+ month_term[month - 1]
+ day;
weekday =
((value * 74899u) >> 16) & 7u;これでも通常のSakamoto式より高速でしたが、月補正値の加算が最終乗算の前にあるため、計算の依存鎖に含まれます。
そこで、乗算の分配法則を使います。
(year term + month term + day) * 74899
= (year term + day) * 74899
+ month term * 74899月補正値は12個しかないため、あらかじめ74899倍した値を保持できます。
通常の月補正値は次のとおりです。
0, 3, 2, 5, 0, 3, 5, 1, 4, 6, 2, 4これを74899倍すると、次の値になります。
0, 224697, 149798, 374495,
0, 224697, 374495, 74899,
299596, 449394, 149798, 299596この変形により、月補正値の加算を最終乗算の後へ移動できます。
Apple Silicon上のClangでは、この部分が積和命令を使ったコードへコンパイルされました。
月テーブルは12バイトから48バイトへ増えますが、年ごとのテーブルや曜日結果のテーブルを追加するわけではありません。一般的なPCやサーバーではキャッシュライン1本に収まる大きさです。メモリやコードサイズを優先する環境では、12バイトの通常月テーブルを使う分配前の実装も選択できます。
なぜ5243で100除算できるのか
補正後の年をyearとします。
対象範囲は次のとおりです。
0 <= year <= 9999yearを次のように分解します。
year = 100 * q + r
0 <= r < 100ここで、
5243 * 100
= 524300
= 2^19 + 12です。
したがって、
5243 * year
= 5243 * (100 * q + r)
= 2^19 * q
+ 12 * q
+ 5243 * rとなります。
year <= 9999なので、q <= 99かつr <= 99です。
残りの部分の最大値は、
12 * 99 + 5243 * 99
= 520245
< 2^19です。
つまり、12 * q + 5243 * rから19ビット目への桁上がりはありません。
そのため、
(year * 5243u) >> 19は、対象範囲で常にyear / 100と一致します。
同じ積で400除算できる理由
今度はyearを次のように分解します。
year = 400 * q + r
0 <= r < 400また、
5243 * 400
= 2097200
= 2^21 + 48です。
したがって、
5243 * year
= 2^21 * q
+ 48 * q
+ 5243 * rとなります。
year <= 9999より、q <= 24かつr <= 399です。
残りの部分の最大値は、
48 * 24 + 5243 * 399
= 2093109
< 2^21です。
したがって、21ビット目への桁上がりはありません。
そのため、
(year * 5243u) >> 21は、対象範囲でyear / 400と一致します。
以上から、一度の乗算結果を共有して次の二つの値を取得できます。
const uint32_t product = year * 5243u;
const uint32_t century = product >> 19;
const uint32_t four_centuries = product >> 21;なぜ74899でmod 7を計算できるのか
Sakamoto式の途中結果をvalueとします。
今回の入力範囲では、valueは次の範囲に収まります。
0 <= value < 13000valueを次のように分解します。
value = 7 * q + r
0 <= r < 7定数74899には、次の関係があります。
7 * 74899
= 8 * 2^16 + 5また、
74899
= 2^16 + 9363です。
これらを使って展開すると、
74899 * value
= 74899 * (7 * q + r)
= (8 * q + r) * 2^16
+ 5 * q
+ 9363 * rとなります。
value < 13000なので、q <= 1857です。
残りの部分の最大値は、
5 * 1857 + 9363 * 6
= 65463
< 2^16です。
したがって、5 * q + 9363 * rから16ビット目への桁上がりはありません。
そのため、
(74899 * value) >> 16
= 8 * q + rとなります。
8 * qの下位3ビットは必ずゼロです。
したがって、下位3ビットを取り出せばrだけが残ります。
((value * 74899u) >> 16) & 7uこれは対象範囲でvalue % 7と一致します。
32-bit整数の範囲
この実装は、32-bit整数のオーバーフローを利用していません。
最初の積の最大値は次のとおりです。
9999 * 5243
= 52424757
< 2^32最終段について、value < 13000という少し広い上限を使っても、
12999 * 74899
= 973612101
< 2^32です。
したがって、すべての積は32-bit符号なし整数に収まります。
全入力での検証
今回のアルゴリズムは年1〜9999に限定しています。
そのため、対象となるすべての入力について、通常のSakamoto式と比較できます。
次はTypeScriptによる全件検証です。
const monthTerm = [0, 3, 2, 5, 0, 3, 5, 1, 4, 6, 2, 4] as const;
function weekdaySakamoto(year: number, month: number, day: number): number {
year -= Number(month < 3);
return (
(year +
Math.floor(year / 4) -
Math.floor(year / 100) +
Math.floor(year / 400) +
monthTerm[month - 1] +
day) %
7
);
}
for (let year = 1; year <= 9_999; year++) {
for (let month = 1; month <= 12; month++) {
for (let day = 1; day <= 31; day++) {
const expected = weekdaySakamoto(year, month, day);
const actual = weekdayBoundedDistributed(year, month, day);
if (actual !== expected) {
throw new Error(`mismatch: ${year}-${month}-${day}: ` + `${actual} !== ${expected}`);
}
}
}
}
console.log("verified");比較する入力数は次のとおりです。
9999 * 12 * 31
= 3,719,628年1〜9999、月1〜12、日1〜31のすべてで、通常のSakamoto式との一致を確認しました。
日については各月の日数を考慮せず、すべて31日まで検証しています。したがって、実在する日付より広い入力集合で一致しています。
この関数自体は、4月31日のような実在しない日付を拒否しません。入力値が実在する日付かどうかの検証は、呼び出し側の責務です。
定数探索
今回採用した式の形について、定数とシフト量も探索しました。
100除算と400除算を同じ積から取得する次の形式では、
product = year * C;
year / 100 = product >> shift;
year / 400 = product >> (shift + 2);年0〜9999の全域で成立する最小シフトは19でした。
その最小シフトにおける定数は、5243のみでした。
また、最終剰余の次の形式では、
((value * C) >> shift) & 7valueの対象範囲全体で成立する最小シフトは16でした。
基本となる定数範囲内では、定数は74899のみでした。
これは、あくまで上記の式の形に限定した最小性です。
あらゆるプログラムやstraight-line expressionの中で、この実装が最適であることを意味するものではありません。
ベンチマーク
Apple Silicon上で、Clangの次の条件を使って測定しました。
-O3 -march=nativeランダムな年、月、日をあらかじめ配列へ格納し、それを繰り返し処理しています。
インライン・バッチ処理の結果は、おおむね次のとおりでした。
| 実装 | 時間 |
|---|---|
| 通常のSakamoto | 約0.528 ns/件 |
| 最初のマジックナンバー版 | 約0.430 ns/件 |
| 月補正分配版 | 約0.370〜0.373 ns/件 |
通常のSakamoto実装に対する月補正分配版の改善率は、約29〜30%でした。
scalar呼び出しでも高速化を確認しましたが、改善率はインライン化の有無、関数呼び出しの方法、コンパイラが生成する命令列によって変わります。
また、Intel x86-64、Rust、Node.js、ブラウザJavaScriptなどでは、同じ改善率になるとは限りません。
掲載したRust版とTypeScript版については、通常のSakamoto式との全件一致を確認しています。性能値を比較する場合は、対象となるCPU、コンパイラ、JITごとに測定する必要があります。
先行研究と今回の位置づけ
5243による0以上9999以下の整数の100除算は既知のstrength reductionです。
例えば、整数から10進文字列への高速変換でも、次の形が使用されています。
(value * 5243) >> 19したがって、定数5243自体を新しいものとして主張するものではありません。
また、2026年8月17日にBen Joffe氏が公開した次の記事では、7が2^3 - 1であることを利用し、限定範囲の剰余を乗算とシフトで求める原理が紹介されています。
今回の実装のポイントは、これらの考え方を年1〜9999に限定したSakamoto式へ適用し、次の一つの実装へまとめたことです。
- 同じ積
year * 5243からyear / 100とyear / 400を取得する - 定数
74899によって最終的なmod 7を3ビット抽出へ変換する - 月補正値を事前にスケーリングする
- 月補正の加算を依存鎖の後段へ移動する
同じ定数構成とSakamoto式への適用については、調査した範囲では先行例を確認できていません。
ただし、構成要素となる定数除算、strength reduction、限定範囲での剰余計算は既知の技法です。
まとめ
年を1〜9999に限定することで、Sakamotoの曜日計算に含まれる定数除算とmod 7を、乗算、シフト、ビット抽出へ置き換えられました。
主な変換は次の二つです。
const uint32_t product = year * 5243u;
product >> 19; // year / 100
product >> 21; // year / 400((value * 74899u) >> 16) & 7uさらに、月補正値をあらかじめ74899倍することで、月項を年計算から続く依存鎖の後段へ移動しました。
手元のApple Silicon環境では、Cのインライン・バッチ処理で通常のSakamoto実装より約30%高速になりました。
一方、この実装には年1〜9999という適用範囲があります。また、性能差はCPUやコンパイラ、JITによって変化します。
一般的なアプリケーションコードでは、読みやすく範囲制限のない通常のSakamoto式で十分です。今回の実装は、曜日計算がホットパスにあり、入力範囲が明確で、性能差を実測できる場合に向いています。