Knuth 4.6.3
fn llpow(r: []Limb, a: []const Limb, b: u32, tmp_limbs: []Limb) void
fn llpow(r: []Limb, a: []const Limb, b: u32, tmp_limbs: []Limb) void {
var tmp1: []Limb = undefined;
var tmp2: []Limb = undefined;
// Multiplication requires no aliasing between the operand and the result
// variable, use the output limbs and another temporary set to overcome this
// limitation.
// The initial assignment makes the result end in `r` so an extra memory
// copy is saved, each 1 flips the index twice so it's only the zeros that
// matter.
const b_leading_zeros = @clz(b);
const exp_zeros = @popCount(~b) - b_leading_zeros;
if (exp_zeros & 1 != 0) {
tmp1 = tmp_limbs;
tmp2 = r;
} else {
tmp1 = r;
tmp2 = tmp_limbs;
}
@memcpy(tmp1[0..a.len], a);
@memset(tmp1[a.len..], 0);
// Scan the exponent as a binary number, from left to right, dropping the
// most significant bit set.
// Square the result if the current bit is zero, square and multiply by a if
// it is one.
const exp_bits = 32 - 1 - b_leading_zeros;
var exp = b << @as(u5, @intCast(1 + b_leading_zeros));
var i: usize = 0;
while (i < exp_bits) : (i += 1) {
// Square
@memset(tmp2, 0);
llsquareBasecase(tmp2, tmp1[0..llnormalize(tmp1)]);
mem.swap([]Limb, &tmp1, &tmp2);
// Multiply by a
const ov = @shlWithOverflow(exp, 1);
exp = ov[0];
if (ov[1] != 0) {
@memset(tmp2, 0);
llmulacc(.add, null, tmp2, tmp1[0..llnormalize(tmp1)], a);
mem.swap([]Limb, &tmp1, &tmp2);
}
}
}