うるう年判定の1073750999はどうやって出てきたのか

Xで、こんなうるう年判定を見かけた。元の投稿は見失ってしまった。（Xの仕様で流れていった）

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

0年から102499年まで閏年判定として使えるらしい。
閏年判定のロジック説明でいうならこれで終わり。

問題は1073750999という数字、はどっからきたの？
年代数えて規則的に出る数字では全然ないので、どうやって導出したのかを調べた方が面白そうなんで、やってみた。

普通の判定式から%25だけを残す

グレゴリオ暦のうるう年は、次の条件。

4で割り切れる
かつ
100で割り切れない、または400で割り切れる

Cで書くとこんな感じ。

bool is_leap_year(uint32_t y) {
    return y % 4 == 0 && (y % 100 != 0 || y % 400 == 0);
}

ここで最初の条件を満たすyは4の倍数。

次に、100は $4 \times 25$ なので、% 100を% 25に置き換えられる。
となると、400は $25 \times 16$ なので、% 16 == 0が% 400 == 0と同じ条件になる。

4と16は2の累乗なので、下位ビットのマスクで剰余を判定できる。

bool target(uint32_t y) {
    return (y & 3) == 0 &&
           (y % 25 != 0 || (y & 15) == 0);
}

これで剰余は% 25だけになる。
こいつがどうやったらああなるのか。

候補式をZ3の制約に合わせる

Z3はSMTソルバーの一つで、整数やビット列の計算まで含めた条件が成立するかを調べられる。
式を成立させる値があればsat、なければunsat、Z3が判定できなければunknownになる。
satなら、条件を満たした具体的な値も取り出せる。

いちいち候補を一つずつ試すコードを書く必要はない。
正解式と候補式が一致する条件を、Z3の制約としてそのまま書ける。
今回なら、その条件からf, m, tを探す。

高速版を考案したFalk Hüffner氏の原典によると、短い式が存在するならマジックナンバーを使ったハッシュのような形になると考えたという。
(y * f) <= tを試した後、マスクを一つ加えた形をまず考える。

((y * f) & m) <= t

そして、空欄に入る乗数f、マスクm、しきい値tをZ3で探す。
乗算1回、AND 1回、比較1回という式の形はHüffner氏がここまで決めてからぶん投げている。

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

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

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

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

Z3には、「32ビット符号なし整数として掛け算をオーバーフローさせ、その結果をマスクしてしきい値と比べる。この式が指定範囲の正解と常に一致するf, m, tは存在するか」という制約を与えている。

ForAllはyを範囲内の全入力として扱う。
Z3がsatを返し、f, m, tの具体値を示した場合、その制約のもとでは正解式と候補式が範囲内の全入力で一致している、ということになる。

18ビット版をZ3で合成した

32ビットの3つの定数を見ても、正直何をしているのか分からない。
なので、まずビット幅を縮めて同じ式を作ってみた。

8ビットでは400年すら表現できない。16ビットで0年から400年までを制約に入れるとunsat。17ビットでやっとsatになった。

18ビットで0年から499年までの500入力を入れたところ、Z3 5.1.0は次の解を返してきた。

f = 65623  = 0x10057
m = 197119 = 0x301FF
t = 496    = 0x001F0

bool leap18(uint32_t y) {
    const uint32_t p = (y * 0x10057u) & 0x3FFFFu;
    return (p & 0x301FFu) <= 0x001F0u;
}

0x3FFFFで18ビット乗算のオーバーフローを再現している。0年から499年までは全件一致し、500年を加えた制約はunsatになった。
この式と18ビット幅では、500年まで正しくなる別の定数は存在しない。

合成と全数照合に使ったコードは LiltingChannelLabo に置いた。

18ビット版の乗算結果を分解する

Z3が返した0x10057を見てても、やっぱり何をしているのか分からない。
まず2項に分けてみる。

0x10057 = 2^16 + 87

18ビット整数で2^16を掛けると、yの下位2ビットが積の最上位2ビットへ移る。

残りのy * 87は、499年までなら $499 \times 87 = 43413 < 2^{16}$ なので、最上位2ビットへ届かない。
積の最上位2ビットが0ならy % 4 == 0になる。

87は、87.04を切り下げた値だった模様。
ここでようやく100が出てきた。

2^9 \times \frac{17}{100} = 87.04

100の倍数を $y=100k$ とすると、積は次の式になる。

y \times 87 = 2^9 \times 17k - 4k

積は $2^9$ の整数倍から $4k$ を引いた値になり、下位側のビットが1で埋まる。

マスクとしきい値を分解すると次の値になる。

m = 0x30000 | 0x001F0 | 0x0000F
t =           0x001F0

0x30000は積の最上位2ビットを取り出し、4の倍数でなければマスク後の値をしきい値より大きくする。
0x001F0は100の倍数が表れる位置の連続5ビットを選ぶ。
0x0000Fは最下位4ビットを使って16の倍数かを判定する。

0年から500年まで、100年刻みの値は次のとおり。

年`y`	`(y * f) mod 2^18`	`& m`	`<= t`	暦の答え
0	`0x00000`	`0x00000`	true	うるう年
100	`0x021FC`	`0x001FC`	false	平年
200	`0x043F8`	`0x001F8`	false	平年
300	`0x065F4`	`0x001F4`	false	平年
400	`0x087F0`	`0x001F0`	true	うるう年
500	`0x0A9EC`	`0x001EC`	true	本当は平年

100、200、300年はマスク後の値がしきい値を少し超える。400年ではちょうどしきい値に重なり、<=なのでうるう年になる。500年になると溜まってた0.04 × yの微妙な数のズレが次のビットへ達して、しきい値より小さくなって誤判定する。

1073750999を16進数へ戻す

ということは同じ手順で32ビット版の1073750999も分解すればわかる。

f = 1073750999 = 0x400023D7 = 2^30 + 9175
m = 3221352463 = 0xC001F00F
t =     126976 = 0x0001F000

Z3が返した9175を100倍すると917500で、 $2^{17}$ の7倍より4小さい。

100 \times 9175 = 7 \times 2^{17} - 4

この式の両辺を100で割ると、次の形になる。

9175 = 2^{17} \times \frac{7}{100} - 0.04

10進数をビット演算の式へ戻すと次の値になる。

f = 2^30 + floor(2^17 * 7 / 100)
m = (3 << 30) | (31 << 12) | 15
t =              (31 << 12)

fの $2^{30}$ は、yの下位2ビットを積の最上位2ビットへ移す。範囲上限でも $102499 \times 9175 = 940428325 < 2^{30}$ なので、9175倍した側から最上位2ビットへの桁上がりはない。

mは、最上位2ビット、中間の連続5ビット、最下位4ビットだけが1になっている。tは中間5ビットだけが1になる。

        最上位2ビット     中間5ビット      最下位4ビット
m = 11 0000000000000 11111 00000000 1111
t = 00 0000000000000 11111 00000000 0000

マスク後の最上位2ビットに1があればtより大きくなり、y % 4 != 0として平年と判定される。
最上位が0で、中間5ビットが全部1でなければtより小さい。これは、4の倍数かつ100の倍数ではない年に対応する。
中間5ビットが全部1なら、最下位4ビットも0のときだけtと等しい。fの下位4ビットは7で奇数なので、積の下位4ビットが0になる条件はy % 16 == 0と等しく、100の倍数のうち400の倍数だけが該当する、という流れ。

9175倍の中間5ビットには、100の倍数以外も混ざる。範囲の前半では下2桁が14、57、71の年も同じ位置に表れるが、4の倍数ではないので最上位2ビットが0にならず、平年と判定される。

9175倍の切り捨て誤差

9175 = 2^17 × 7 / 100 - 0.04なので、年が増えると切り捨てた0.04の誤差も増える。102500年では、中間5ビットが100の倍数を表す位置からズレた。

年`y`	`(y * f) mod 2^32`	`& m`	`<= t`	暦の答え
102499	`0xF80DCC25`	`0xC001C005`	false	平年
102500	`0x380DEFFC`	`0x0001E00C`	true	本当は平年

102500年は100で割り切れて400では割り切れない平年だが、マスク後の値がしきい値より小さくなり、うるう年と誤判定する。
全数照合で確認してもこれが最初の不一致で、102499という数字はその直前の値だった。