I read A faster way to convert a timestamp ➜ Hour, Min, Sec the other day and coincidentally I had been messing around with another problem involving the multiply by a magic number instead of division trick. In my problem I was trying to do this for a u16 in SIMD. AVX2 has mulhi_epu16 but not mulhi_epu32 so then you have to go all the way to mul_epu32 which takes us from 16 lanes down to just 4 if we have to do a full 32x32->64 bit multiply.

So I was tinkering around trying to see what the smallest magic number you can use is. You can read the LLVM source for division by a constant, but it didn’t answer my question, so I reached for the Z3 hammer.

The problem I was playing around with was actually inspired after reading about the Go Greentea garbage collector, and then seeing in the code that they use a magic number with no shift, whereas in general you need a shift. See the source here. This part wasn’t added by Greentea but it had my mind on SIMD. Go uses this magic shift to turn a u16 page offset into a slot index by dividing by the u16 size class. I wanted to do this in SIMD hence the intro of the article. Turns out, you need a mulhi32 to do this with no shift. You can do it with mulhi16 and a single shift per size class. This is slightly less ideal in AVX2 because there is no srlv_epi16, but could be done with srlv_epi32 and masking. Note that turning offsets into slots is maybe not that useful because once you have slot indices, you then have a gather problem which isn’t good with AVX2.

Here is the Z3 code that I experimented with (note that in Z3 ZeroExt(n, x) adds n zeros, not bringing it to n bits):

from z3 import Solver, BitVec, LShR, UDiv, BitVecVal, ForAll, Extract, ZeroExt, sat, Not, And, Exists
import math

def hi(x):
    n = x.size()
    return Extract(n-1, n//2, x)

def lo(x):
    n = x.size()
    return Extract(n//2 - 1, 0, x)

def mulhi(a, b):
    assert a.size() == b.size(), (a.size(), b.size())
    n = a.size()
    return hi(ZeroExt(n, a) * ZeroExt(n, b))

def ZeroExtTo(n, x):
    assert x.size() <= n
    if x.size() < n:
        return ZeroExt(n - x.size(), x)
    return x

def format_general_model(d, m, s1, s2):
    if s1 == 0 and s2 == 0:
        return f'mulhi(a, {m}) == a // {d}'
    if s1 == 0:
        return f'mulhi(a, {m}) >> {s2} == a // {d}'
    if s2 == 0:
        return f'mulhi(a >> {s1}, {m}) == a // {d}'
    return f'mulhi(a >> {s1}, {m}) >> {s2} == a // {d}'

# realizing this doesn't include possible adds
def solve_general(Bx, Bm, divisor, *, Bd=None, no_s1=False, no_s2=False):
    if Bd is None:
        Bd = Bx
    assert divisor < 2**Bd
    B = max(Bx, Bm, Bd)
    if B % 2 == 1:
        B += 1
    Bs = int(math.ceil(math.log2(B)))
    s = Solver()
    x = BitVec('x', Bx)
    d = BitVecVal(divisor, Bd)
    m = BitVec('m', Bm)
    s1 = BitVec('s1', Bs)
    s2 = BitVec('s2', Bs)

    zx = ZeroExtTo(B, x)
    zd = ZeroExtTo(B, d)
    zm = ZeroExtTo(B, m)
    zs1 = ZeroExtTo(B, s1)
    zs2 = ZeroExtTo(B, s2)

    # mulhi(x >> s1, m) >> s2 == x // d
    s.add(ForAll([x], LShR(mulhi(LShR(zx, zs1), zm), zs2) == UDiv(zx, zd)))

    if no_s1:
        s.add(s1 == 0)
    if no_s2:
        s.add(s2 == 0)

    while s.check() == sat:
        model = s.model()
        yield model[m], model[s1], model[s2]

        s.add(Not(And(m == model[m], s1 == model[s1], s2 == model[s2])))

go_size_classes = [8, 16, 24, 32, 48, 64, 80, 96, 112, 128, 144, 160, 176, 192, 208, 224, 240, 256, 288, 320, 352, 384, 416, 448, 480, 512, 576, 640, 704, 768, 896, 1024, 1152, 1280, 1408, 1536, 1792, 2048, 2304, 2688, 3072, 3200, 3456, 4096, 4864, 5376, 6144, 6528, 6784, 6912, 8192, 9472, 9728, 10240, 10880, 12288, 13568, 14336, 16384, 18432, 19072, 20480, 21760, 24576, 27264, 28672, 32768]

Bx = 16

print('-- Bm == 16 --')
for size in go_size_classes:
    for model in solve_general(Bx, Bx, size):
        print(format_general_model(size, *model))
        break
    else:
        print(size, 'unsat')

print('\n-- Bm == 16 and no s1 --')
for size in go_size_classes:
    for model in solve_general(Bx, Bx, size, no_s1=True):
        print(format_general_model(size, *model))
        break
    else:
        print(size, 'unsat')

print('\n-- Bm == 16 and no s2 --')
for size in go_size_classes:
    for model in solve_general(Bx, Bx, size, no_s2=True):
        print(format_general_model(size, *model))
        break
    else:
        print(size, 'unsat')

print('\n-- Bm == 16 and no s1 or s2 --')
for size in go_size_classes:
    for model in solve_general(Bx, Bx, size, no_s1=True, no_s2=True):
        print(format_general_model(size, *model))
        break
    else:
        print(size, 'unsat')

print('\n-- Bm == 32 and no shift --')
for size in go_size_classes:
    for model in solve_general(Bx, Bx*2, size, no_s1=True, no_s2=True):
        print(format_general_model(size, *model))
        break
    else:
        print('unsat')

One interesting thing to note is that the shifts and magic numbers aren’t necessarily unique.

Anyways, looping back to the timestamp -> hh:mm:ss conversion: you’re given a timestamp in number of seconds and want to turn it into a clock timestamp. The linked article shows an algorithm that does two division with two multiplies right away, dividing t / 60 and t / 3600 and then some followup calculations to get your result. I then wondered if you could do both divides (via multiplies) in a SWAR fashion.

Because Z3 was fresh in my head, I wrote this script to find the magic number and shifts for just that:

from z3 import Solver, BitVec, LShR, UDiv, BitVecVal, ForAll, Extract, ZeroExt, sat, Not, And

def hi(x):
    n = x.size()
    return Extract(n-1, n//2, x)

def lo(x):
    n = x.size()
    return Extract(n//2 - 1, 0, x)

Bx = 16
b = 6

d1 = 60
d2 = 3600

s = Solver()
x = BitVec('x', Bx)
m = BitVec('m', 64)
s1 = BitVec('s1', b)
s2 = BitVec('s2', b)
s3 = BitVec('s3', b)
s4 = BitVec('s4', b)

zx = ZeroExt(128 - Bx, x)
zm = ZeroExt(64, m)
zs1 = ZeroExt(64-b, s1)
zs2 = ZeroExt(64-b, s2)
zs3 = ZeroExt(64-b, s3)
zs4 = ZeroExt(64-b, s4)

q1 = ZeroExt(64 - Bx, UDiv(x, d1))
q2 = ZeroExt(64 - Bx, UDiv(x, d2))

zxm = zx * zm

ones = BitVecVal(-1, 64)

q1_ = LShR(lo(zxm), zs1) & LShR(ones, zs3)
q2_ = LShR(hi(zxm), zs2) & LShR(ones, zs4)

s.add(ForAll([x], And(q1 == q1_, q2 == q2_)))

min_s3 = 64
min_s4 = 64
s.push()
s.add(And(s3 < min_s3, s4 < min_s4))

while s.check() == sat:
    model = s.model()
    m_ = model[m].as_long()
    s1_ = model[s1].as_long()
    s2_ = model[s2].as_long()
    s3_ = model[s3].as_long()
    s4_ = model[s3].as_long()
    print(f'0x{m_:016x} s1={s1_} s2={s2_} s3={s3_} ({64 - s3_} bits) s4={s4_} ({64 - s4_} bits)')

    min_s3 = min(min_s3, s3_)
    min_s4 = min(min_s4, s4_)
    s.pop()
    s.push()
    s.add(And(s3 < min_s3, s4 < min_s4))

This computes a 64 bit multiplier m s.t. (mullo(m, x) >> s1) & (-1 >> s3) == x // d1 and (mulhi(m, x) >> s2) & (-1 >> s4). The mask with -1 >> s says the low (64 - s) bits are correct. Ideally we want this to be the most amount of bits, so those shifts should be as small as possible. Here I just iteratively lower the bound until we get unsat. Maybe some scenarios you want one answer with more bits than the other (foreshadowing) and you could find other solutions. For example, the above prints out:

0x091a347000022223 s1=23 s2=7 s3=53 (11 bits) s4=53 (11 bits)
0x48d159f000044445 s1=24 s2=10 s3=52 (12 bits) s4=52 (12 bits)
0x1234685800008889 s1=21 s2=8 s3=51 (13 bits) s4=51 (13 bits)
0x2468c00000044446 s1=24 s2=9 s3=50 (14 bits) s4=50 (14 bits)
0x91a2e00000444465 s1=28 s2=11 s3=47 (17 bits) s4=47 (17 bits)
0x48d15a0000008889 s1=21 s2=10 s3=46 (18 bits) s4=46 (18 bits)
0x48d1a00000022223 s1=23 s2=10 s3=45 (19 bits) s4=45 (19 bits)
0x48d1800000044446 s1=24 s2=10 s3=41 (23 bits) s4=41 (23 bits)
0x91a3000000044447 s1=24 s2=11 s3=40 (24 bits) s4=40 (24 bits)
0x91a3400000008889 s1=21 s2=11 s3=39 (25 bits) s4=39 (25 bits)
0x48d1800000008889 s1=21 s2=10 s3=38 (26 bits) s4=38 (26 bits)
0x91a3000000008889 s1=21 s2=11 s3=37 (27 bits) s4=37 (27 bits)

Another consideration with how many bits we get is that 32 bits (or at least 32 bits) of correctness is very convenient because we don’t need a mask (on x86-64).

In writing this post, I discovered why something I had tried earlier worked accidentally and is a nice bonus. In doing our conversion, we compute:

to compute:
hh = t / 3600
mm = (t % 3600) / 60
ss = t % 60

identity 1:
x % D = (x + c*(x/D)) % (D + c)

rewrite :
hh = t / 3600

mm = (t % 3600) / 60
   = (t/60) % 60
   = ((t/60) + 4 * ((t/60) / 60)) % (60 + 4)
   = ((t/60) + 4 * (t/3600)) % 64

ss = t % 60
   = (t + 4 * (t/60)) % 64

compute:
div3600 = t / 3600
div60 = t / 60

hh = div3600
mm = (div60 + 4 * div3600) & 63
ss = (t + 4 * div60) & 63

Notice that in both uses of div60 we are taking the answer mod 64, so we actually only care about div60 being correct to 6 bits! That is great because it means we don’t have to mask after we shift. It also means we can find a magic number that gives us more correct bits for div3600.

from z3 import Solver, BitVec, LShR, UDiv, BitVecVal, ForAll, Extract, ZeroExt, sat, Not, And, URem

def hi(x):
    n = x.size()
    return Extract(n-1, n//2, x)

def lo(x):
    n = x.size()
    return Extract(n//2 - 1, 0, x)

Bx = 20
b = 6

s = Solver()
x = BitVec('x', Bx)
m = BitVec('m', 64)
s1 = BitVec('s1', b)
s2 = BitVec('s2', b)
s3 = BitVec('s3', b)
s4 = BitVec('s4', b)
zx = ZeroExt(128 - Bx, x)
zm = ZeroExt(64, m)
zs1 = ZeroExt(64 - b, s1)
zs2 = ZeroExt(64 - b, s2)
zs3 = ZeroExt(64 - b, s3)

ones = BitVecVal(-1, 64)

q1 = ZeroExt(64 - Bx, URem(UDiv(x, 60), 64))
q2 = ZeroExt(64 - Bx, UDiv(x, 3600))

zxm = zx * zm

q1_ = LShR(lo(zxm), zs1) & BitVecVal(63, 64)
q2_ = LShR(hi(zxm), zs2) & LShR(ones, zs3)

s.add(ForAll([x], And(q1 == q1_, q2 == q2_)))

min_s3 = 33

s.push()
s.add(s3 < min_s3)

while s.check() == sat:
    model = s.model()
    m_ = model[m].as_long()
    s1_ = model[s1].as_long()
    s2_ = model[s2].as_long()
    s3_ = model[s3].as_long()
    print(f'0x{m_:016x} s1={s1_} s2={s2_} s3={s3_} q2 correct to {64 - s3_} bits')
    min_s3 = min(min_s3, s3_)
    s.pop()
    s.push()
    s.add(s3 < min_s3)

After fiddling around with Bx the number of input bits we are correct to (1 full day is 16.4 bits), I happily discovered a solution at Bx=20 where s2=0. This means we don’t even have to shift that half! Another solution is shown that has a slightly higher range. Higher values of Bx are maybe possible but Z3 gets a bit slow and I didn’t wait around. (To be honest I’m not sure the practicality of this problem statement with either t in 1 day already or not also taking hh % 24)

struct TIME_S {
    uint32_t hour, min, sec;
};

// Accurate to 1198378
TIME_S time_bogo_1(uint32_t t) {
    uint128_t prd = (uint128_t)t * 0x0012345700111112;
    uint64_t lo = prd;
    uint64_t hi = prd >> 64;
    // here div3600 is in the hi bits and div60 is in the lo bits
    uint32_t div60 = lo >> 26;

    uint32_t hh = hi; // no shift!
    uint32_t mm = (div60 + 4 * hh) % 64;
    uint32_t ss = (t + 4 * div60) % 64;

    return { hh, mm, ss };
}

// Accurate to 6100858
TIME_S time_bogo_2(uint32_t t) {
    uint128_t prd = (uint128_t)t * 0x91a2b49800444445;
    uint64_t lo = prd;
    uint64_t hi = prd >> 64;
    // here div3600 is in the hi bits and div60 is in the lo bits
    uint32_t div60 = lo >> 28;

    uint32_t hh = hi >> 11;
    uint32_t mm = (div60 + 4 * hh) % 64;
    uint32_t ss = (t + 4 * div60) % 64;

    return { hh, mm, ss };
}

Benchmarking these functions against benjoffe/fast-world-calendars on a 5950x gives (many rows deleted):

> build/bin/time_bench
----------------------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations UserCounters...
----------------------------------------------------------------------------------
bench_baseline               19851 ns        19730 ns        35481 ns/t=0.3029
bench_bogo_1                 62663 ns        62227 ns        11254 ns/t=0.956172
bench_bogo_2                 68966 ns        68508 ns        10218 ns/t=1.05234
bench_benjoffe_v1            78708 ns        78167 ns         8957 ns/t=1.20099
bench_benjoffe_v2           107858 ns       107152 ns         6527 ns/t=1.64579
bench_benjoffe_v3            78545 ns        77997 ns         8974 ns/t=1.1985
bench_traditional1          117706 ns       116864 ns         5991 ns/t=1.79605
bench_traditional1_fp        81920 ns        81219 ns         8687 ns/t=1.25

> build/bin/time_bench -latency
----------------------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations UserCounters...
----------------------------------------------------------------------------------
bench_baseline                9852 ns         9798 ns        71843 ns/t=0.150326
bench_bogo_1                176188 ns       175081 ns         3998 ns/t=2.68843
bench_bogo_2                195286 ns       194348 ns         3602 ns/t=2.97984
bench_benjoffe_v3           195393 ns       194398 ns         3601 ns/t=2.98147
bench_benjoffe_v3b          195357 ns       194344 ns         3601 ns/t=2.98092
bench_benjoffe_v3_64        195391 ns       194331 ns         3602 ns/t=2.98145
bench_benjoffe_v3_64b       195402 ns       194348 ns         3602 ns/t=2.98162

The range isn’t nearly as high as other methods, but I think doing it in one multiply is pretty fun.

The asm for time_bogo_1 (godbolt):

time_bogo_1(unsigned int):
        movabs  rax, 5124097848709394
        mov     edx, edi
        mulx    rcx, rdx, rax
        shr     rdx, 26
        lea     eax, [rdx + 4*rcx]
        lea     edx, [rdi + 4*rdx]
        and     eax, 63
        and     edx, 63
        shl     rax, 32
        or      rax, rcx ; returns in RAX,RDX
        ret

In case you’re not familiar, BOGO is “buy-one-get-one”.