How Z3 found the leap-year magic number 1073750999

I saw this leap-year test on X, then lost the original post. It disappeared into the X timeline.

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

Apparently, it works as a leap-year test from year 0 through 102499.
If all I wanted to explain was the test itself, that would be the end.

The problem is the number 1073750999. Where did that come from?
Counting through the years is never going to produce a number that looks like this, so figuring out how it was derived sounded more interesting. I tried.

Reduce the ordinary test to %25

The Gregorian leap-year rule is:

divisible by 4
and
not divisible by 100, or divisible by 400

In C, that looks like this:

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

At this point, y is already a multiple of 4.

Since $100 = 4 \times 25$, % 100 can be replaced with % 25.
Then $400 = 25 \times 16$, so % 16 == 0 is equivalent to % 400 == 0 here.

Both 4 and 16 are powers of two, so those remainders can be tested with masks over the low bits.

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

Now % 25 is the only remainder operation left.
How does that turn into the strange expression at the top?

Turn the candidate expression into a Z3 constraint

Z3 is an SMT solver. It can decide whether a condition involving integer or bit-vector operations can hold.
If values satisfying the formula exist, the result is sat; if none exist, it is unsat; and if Z3 cannot decide, it returns unknown.
For a sat result, the concrete values that satisfy the condition can also be retrieved.

There is no need to write a loop that tries candidate constants one at a time.
The equality between the correct expression and a candidate expression can be written directly as a Z3 constraint.
Here, the unknown values are f, m, t.

According to Falk Hüffner’s original article, he expected any short expression to look like a weird hash built from magic constants.
After trying (y * f) <= t, he added one mask and started with this form:

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

Z3 then searches for the multiplier f, mask m, and threshold t.
Hüffner had already decided that the expression would use one multiply, one AND, and one comparison before throwing the remaining work at Z3.

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)),
))

The constraint says: multiply as an unsigned 32-bit integer and allow it to overflow, mask the result, and compare it with a threshold. Do any values of f, m, t make that expression agree with the correct answer throughout the requested range?

ForAll treats y as every input in the range.
If Z3 returns sat and supplies concrete values for f, m, t, the correct and candidate expressions agree for every input covered by that constraint.

Synthesize an 18-bit version with Z3

Three huge 32-bit constants do not tell me much about what the expression is doing.
So I reduced the bit width and built the same expression first.

Eight bits cannot even represent a full 400-year cycle. With constraints for years 0 through 400, 16 bits gives unsat, while 17 bits finally gives sat.

For an 18-bit version constrained over the 500 inputs from year 0 through 499, Z3 5.1.0 returned:

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;
}

The 0x3FFFF mask reproduces overflow in an 18-bit multiplication. The result agrees for every year from 0 through 499, while adding year 500 makes the constraint unsat.
With this expression and an 18-bit width, no set of constants can remain correct through year 500.

The synthesis and exhaustive-checking code is in LiltingChannelLabo.

Split the 18-bit multiplication apart

Looking at the 0x10057 returned by Z3 still did not tell me what it was doing.
I split it into two terms.

0x10057 = 2^16 + 87

Multiplying by 2^16 in an 18-bit integer moves the lowest two bits of y into the highest two bits of the product.

For years up to 499, the remaining y * 87 term satisfies $499 \times 87 = 43413 < 2^{16}$, so it cannot reach those highest two bits.
If the highest two bits of the product are zero, then y % 4 == 0.

The 87 appears to be 87.04 rounded down.
This is where 100 finally shows up.

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

For a multiple of 100, write $y=100k$. The product becomes:

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

The product is a multiple of $2^9$ minus $4k$. That subtraction fills the lower-side bits with ones.

The mask and threshold split into these fields:

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

0x30000 extracts the highest two bits of the product. If the year is not divisible by 4, it makes the masked value greater than the threshold.
0x001F0 selects the five consecutive bits where multiples of 100 appear.
0x0000F uses the lowest four bits to test divisibility by 16.

Here are the values at 100-year intervals from year 0 through 500:

Year y	(y * f) mod 2^18	& m	<= t	Calendar answer
0	0x00000	0x00000	true	leap year
100	0x021FC	0x001FC	false	common year
200	0x043F8	0x001F8	false	common year
300	0x065F4	0x001F4	false	common year
400	0x087F0	0x001F0	true	leap year
500	0x0A9EC	0x001EC	true	actually a common year

At years 100, 200, and 300, the masked value sits just above the threshold. At year 400 it lands exactly on the threshold, so <= is true and the result is a leap year. At year 500, the accumulated 0.04 × y drift has reached the next bit. The masked value drops below the threshold and produces the wrong answer.

Convert 1073750999 back to hexadecimal

The 32-bit 1073750999 should make sense after splitting it the same way.

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

Multiplying Z3’s 9175 by 100 gives 917500, which is four less than seven times $2^{17}$.

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

Dividing both sides by 100 gives the form containing 7.

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

Written as bit operations, the decimal constants become:

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

The $2^{30}$ term in f moves the lowest two bits of y into the highest two bits of the product. Even at the top of the valid range, $102499 \times 9175 = 940428325 < 2^{30}$, so the 9175 * y side never carries into those highest two bits.

In m, only the highest two bits, a consecutive five-bit field in the middle, and the lowest four bits are set. In t, only the middle five-bit field is set.

        highest 2 bits     middle 5 bits      lowest 4 bits
m = 11 0000000000000 11111 00000000 1111
t = 00 0000000000000 11111 00000000 0000

If either of the highest two bits survives the mask, the result is greater than t. The test classifies y % 4 != 0 as a common year.
If the highest field is zero and the middle five bits are not all ones, the result is less than t. This corresponds to a multiple of 4 that is not a multiple of 100.
If all five middle bits are one, the result equals t only when the lowest four bits are also zero. The lowest four bits of f are 7, which is odd, so the product is zero modulo 16 exactly when y % 16 == 0. Among multiples of 100, that selects only multiples of 400.

The middle five bits of the 9175 * y product also match some years that are not multiples of 100. Early in the range, years ending in 14, 57, and 71 appear there too. None is divisible by 4, so the highest two bits reject them as common years.

The truncation error in the 9175 term

Because 9175 = 2^17 × 7 / 100 - 0.04, the truncated 0.04 error accumulates as the year grows. At year 102500, the middle five bits move away from the position that represents multiples of 100.

Year y	(y * f) mod 2^32	& m	<= t	Calendar answer
102499	0xF80DCC25	0xC001C005	false	common year
102500	0x380DEFFC	0x0001E00C	true	actually a common year

Year 102500 is divisible by 100 but not by 400, so it is a common year. The masked value falls below the threshold and the expression misclassifies it as a leap year.
An exhaustive check finds this as the first mismatch. Year 102499 is immediately before it.