64-bit Solinas prime Fp modular multiplication

solinas_primes.txt

Reduction as in "Generalized Mersenne Numbers" J. Solinas (NSA).

Prime from "Solinas primes of small weight for fixed sizes" J. Angel.

Fp = 2^64-2^32+1

t = 2^32
d = 2
f(t) = t^2 - c_1 * t - c_2

c_1 = 1
c_2 = -1

X_0j = c_(d-j)
X_00 = c_2
x_01 = c_1

X_i0 = X_(i-1,d-1) * c_d
X_ij = X_(i-1,j-1) + X_(i-1,d-1) * c_(d-j)

X_00 = -1
x_01 = 1

X_i0 = -X_(i-1,1)
X_i1 = X_(i-1,0) + X_(i-1,1)

X_10 = -1
X_11 = 0

[ t^2 t^3 ] = [ [-1 1] [-1 0] ] * [1 t] (mod f(t))

[A0 A1] [1 t] + [A2 A3] [t^2 t^3]
= ([A0 A1] + [A2 A3] * [[-1 1][-1 0]]) * [1 t]

B0 = A0 - A2 + A3
B1 = A1      + A2

R = T + S - D

2^k = 2^32

T = A0 + A1 * 2^k
S = A3 + A2 * 2^k
D = A2

Reduction:

#define ROL64(x, r) ((x >> (32 - (r))) | (x << (r)))

p = (uint128_t)x * y;
p0 = (uint64_t)p;
p1 = (uint64_t)(p >> 64);
uint64_t temp = (p1 >> 32) + (p1 << 32) - (uint32_t)p1;
r = p0 + temp;
if (r < temp) r += (1ull << 32) + 1;