Function sqrt [src]

r = ⌊√a⌋ r may alias a. Asserts that r has enough limbs to store the result. Upper bound is (a.limbs.len - 1) / 2 + 1. limbs_buffer is used for temporary storage. The amount required is given by calcSqrtLimbsBufferLen.

Prototype

pub fn sqrt( r: *Mutable, a: Const, limbs_buffer: []Limb, ) void

Parameters

r: *Mutablea: Constlimbs_buffer: []Limb

Source

  pub fn sqrt(
    r: *Mutable,
    a: Const,
    limbs_buffer: []Limb,
) void {
    // Brent and Zimmermann, Modern Computer Arithmetic, Algorithm 1.13 SqrtInt
    // https://members.loria.fr/PZimmermann/mca/pub226.html
    var buf_index: usize = 0;
    var t = b: {
        const start = buf_index;
        buf_index += a.limbs.len;
        break :b Mutable.init(limbs_buffer[start..buf_index], 0);
    };
    var u = b: {
        const start = buf_index;
        const shift = (a.bitCountAbs() + 1) / 2;
        buf_index += 1 + ((shift / limb_bits) + 1);
        var m = Mutable.init(limbs_buffer[start..buf_index], 1);
        m.shiftLeft(m.toConst(), shift); // u must be >= ⌊√a⌋, and should be as small as possible for efficiency
        break :b m;
    };
    var s = b: {
        const start = buf_index;
        buf_index += u.limbs.len;
        break :b u.toConst().toMutable(limbs_buffer[start..buf_index]);
    };
    var rem = b: {
        const start = buf_index;
        buf_index += s.limbs.len;
        break :b Mutable.init(limbs_buffer[start..buf_index], 0);
    };

    while (true) {
        t.divFloor(&rem, a, s.toConst(), limbs_buffer[buf_index..]);
        t.add(t.toConst(), s.toConst());
        u.shiftRight(t.toConst(), 1);

        if (u.toConst().order(s.toConst()).compare(.gte)) {
            r.copy(s.toConst());
            return;
        }

        // Avoid copying u to s by swapping u and s
        const tmp_s = s;
        s = u;
        u = tmp_s;
    }
}  