Fast 64-bit almost-finite field

fast_field_65.cpp

#include <stdint.h>

/*
    65-bit multiplication of two 64-bit values a, b
    with low bits set to 1 (implicitly):

    (2a + 1)(2b + 1)
    = 4ab + 2a + 2b + 1

    Stuffing it back into a 64-bit value:
    = 2ab + a + b
*/

// Product of a and b
static uint64_t mul65(uint64_t a, uint64_t b)
{
    return ((a * b) << 1) + (a + b);
}

// Multiplicative inverse of a
static uint64_t inv65(uint64_t a)
{
    uint64_t x = mul65(1, a) ^ 1; // 5 bit accurate guess
    x = mul65(x, 1 + ~mul65(a, x)); // 10 bits
    x = mul65(x, 1 + ~mul65(a, x)); // 20 bits
    x = mul65(x, 1 + ~mul65(a, x)); // 40 bits
    x = mul65(x, 1 + ~mul65(a, x)); // 80 bits
    return x;
}

// Sum of a, b, c
static uint64_t add65(uint64_t a, uint64_t b, uint64_t c)
{
    return a + b + c + 1;
}

// = a + b - c
static uint64_t sub65(uint64_t a, uint64_t b, uint64_t c)
{
    return a + b - c;
}