Mnożenie dwóch liczb 32-bitowych za pomocą rejestrów tylko 16-bitowych

Próbuję pomnożyć dwie liczby 32-bitowe, ale używając do tego intów tylko 16-bitowych:

#include <cstdint>
#include <cstdio>
#include <bitset>
#include<iostream>

uint16_t multiply_hi(uint16_t a, uint16_t b)
{
    uint16_t a_lo = (uint8_t)a;
    uint16_t a_hi = a >> 8;
    uint16_t b_lo = (uint8_t)b;
    uint16_t b_hi = b >> 8;

    uint16_t a_x_b_hi = a_hi * b_hi;
    uint16_t a_x_b_mid = a_hi * b_lo;
    uint16_t b_x_a_mid = b_hi * a_lo;
    uint16_t a_x_b_lo = a_lo * b_lo;

    uint16_t carry_bit = ((uint16_t)(uint8_t)a_x_b_mid + (uint16_t)(uint8_t)b_x_a_mid + (a_x_b_lo >> 8)) >> 8;

    uint16_t multhi = a_x_b_hi + (a_x_b_mid >> 8) + (b_x_a_mid >> 8) + carry_bit;

    return multhi;
}

uint16_t add_hi(uint16_t a, uint16_t b)
{
    uint16_t a_lo = (uint8_t)a;
    uint16_t a_hi = a >> 8;
    uint16_t b_lo = (uint8_t)b;
    uint16_t b_hi = b >> 8;

    uint16_t a_b_carry_low = (a_lo + b_lo) >> 8;
    uint16_t a_b_carry_high = (a_hi + b_hi) >> 8;

    uint16_t addhi = (a_b_carry_low + a_b_carry_high);
    if (addhi > 1) { addhi = 1; }

    return addhi;
}

int main()
{
    uint16_t a_low = (uint16_t)4154247521;
    uint16_t a_hi = 4154247521 >> 16;
    uint16_t b_low = (uint16_t)3030958429;
    uint16_t b_hi = 3030958429 >> 16;

    uint16_t wynik_hi = multiply_hi(a_low, b_low) * a_hi * b_hi;
    uint16_t wynik_low = a_low * b_low;

    uint32_t wynik = (wynik_hi << 16) | wynik_low;

    std::cout << wynik << "\n";
}

Funkcja multiply_hi prawidłowo oblicza górne 16 bitów, będące wynikiem mnożenia dwóch 16-bitowych liczb. Ale ja mnożę 32-bitowe liczby, więc to nie jest kompletny wynik. Wydawało mi się, że trzeba to pomnożyć jeszcze przez górne bity obu 32-bitowych liczb, żeby otrzymać ostatecznie 16 górnych bitów wyniku. Ale wynik_hi jest błędny.

Co robię źle? Co jest nie tak z obliczaniem górnych bitów wyniku?

Docelowo chcę mnożyć liczby 256-bitowe używając intów 128-bitowych. Piszę własne funkcje bo np. biblioteka boost multiprecision w przypadku intów 256-bitowych jest 400 razy wolniejsza, niż zwykłe mnożenie intów 64-bitowych, jest absurdalnie wolna.

ale przecież wynik_hi to powinno być (x_lo*y_lo)>>16 + x_lo*y_hi + x_hi*y_lo...