Karatsuba Square Root implementation correctness

Question

I'm currently trying to implement Karatsuba Sqrt for my BigInteger module in luau, though I'm having trouble with its accuracy and I don't know what I'm doing wrong.

Karatsuba function:

local function KaratsubaSqrt(n: AptInt): (AptInt, AptInt)
    --> NOTE: even tho we dont have burnikel-ziegler division, this is still faster by a ton

    -- https://gmplib.org/manual/Square-Root-Algorithm
    -- https://en.wikipedia.org/wiki/Integer_square_root#Karatsuba_square_root_algorithm
        -- https://www.isa-afp.org/browser_info/current/AFP/Karatsuba_Sqrt/outline.pdf

        -- must be positive integer
    if n.signum ~= 1 then
        return AptInt.new(), AptInt.new() -- return 0, 0
    end
    
    --> needed here too, fallback to base case
    local limbLen: number = #n.limbs
    if limbLen < SQRT_KARATSUBA_THRESHOLD then
        local sq = n:sqrt()
        return sq, n:SubtractRaw(sq:MultiplyRaw(sq))
    end
    
    local AP_ONE: AptInt = AptInt.new({1})

    --> calculate a3, a2, a1 and a0
    local m2: number = limbLen // 2
    local m: number = m2 // 2
    
    local a3a2: AptInt, a1a0: AptInt = GetKaratsubaUpper(n, m2), GetKaratsubaLower(n, m2)
    local a1: AptInt, a0: AptInt = GetKaratsubaUpper(a1a0, m), GetKaratsubaLower(a1a0, m)
    
    local s1: AptInt, r1: AptInt = KaratsubaSqrt(a3a2)
        
    --> divrem (r1*b + a1, 2*s1)
    local q, u: AptInt = 
        r1:LeftShift(m):AddRaw(a1)
        :DivideRaw(s1:AddRaw(s1))

    --> reconstruct sqrt
    local s: AptInt = s1:LeftShift(m):AddRaw(q) -- s=s1*b+q

    local r: AptInt = n:SubtractRaw(s:MultiplyRaw(s)) -- or u*b+a0-q^2
    
        --> correction
    if r.signum == -1 then
        r = r:AddRaw(s:AddRaw(s)):SubtractRaw(AP_ONE)
        s = s:SubtractRaw(AP_ONE)
    end

    return s, r
end

Base case sqrt:

function AptInt:sqrt(): AptInt
    --> we compute the isqrt via Newton-Heron iteration

    if self.signum ~= 1 then
        return AptInt.new()
    end
    

    local prev2: AptInt = AP_NEG_ONE:clone()
    local prev1: AptInt = AP_ONE:LeftShift(#self.limbs//2+1)

    while true do
        local x1: AptInt = prev1:AddRaw(self:DivideRaw(prev1)):DivideRaw(AP_TWO)

        local equals: boolean = x1:EqualsRaw(prev1)
        if equals then
            return x1
        end

        if x1:EqualsRaw(prev2) and not equals then
            return (x1:LowerThanRaw(prev1) and x1 or prev1) -- min(x1, prev1)
        end

        prev2, prev1 = prev1, x1
    end
end

reproducible example:

local AptInt = {}
AptInt.__index = AptInt

export type ValidConstructor = string | number | AptInt | {number} | nil

export type AptInt = setmetatable<{
    limbs: {number},
    signum: number,

    read clone: (n: AptInt) -> AptInt,

    read Negate: (n: AptInt, inPlace: boolean) -> AptInt,

    read LeftShift: (n: AptInt, amount: number) -> AptInt,

    read AddRaw: (term: AptInt, term: AptInt) -> AptInt,
    read SubtractRaw: (term: AptInt, term: AptInt) -> AptInt,
    read MultiplyRaw: (factor: AptInt, factor: AptInt) -> AptInt,
    read DivideRaw: (dividend: AptInt, divisor: AptInt) -> (AptInt, AptInt),
    read sqrt: (n: AptInt) -> AptInt,

    read EqualsRaw: (self: AptInt, num: AptInt) -> boolean,
    read LowerThanRaw: (self: AptInt, num: AptInt) -> boolean,
    read LowerOrEqualToRaw: (self: AptInt, num: AptInt) -> boolean,

    read ToString: (n: AptInt) -> string,
}, typeof(AptInt)>

-- ================================
--> OPTIMIAZTION FUNCTION CONSTANTS
-- ================================

local log10: (number) -> number = math.log10
local log: (number, number?) -> number = math.log
local floor: (number) -> number = math.floor
local ceil: (number) -> number = math.ceil
local max: (number) -> number = math.max
local abs: (number) -> number = math.abs
local sign: (number) -> number = math.sign

local tcreate: <V>(number, V?) -> {V} = table.create
local tinsert: <V>({V}, V) -> () & <V>({V}, number, V) -> () = table.insert
local tremove: <V>({V}, number?) -> V? = table.remove
local tclone: <T>({T}) -> {T} = table.clone

local sformat: (string, ...unknown) -> string = string.format

-- =================
--> NORMAL CONSTANTS
-- =================

local SQRT_KARATSUBA_THRESHOLD: number = 16
local SAFE_UNPACK_THRESHOLD: number = 8000 -- unpack() max size in limbs

local BASE_POW: number = 7 -- MUST SATISFY: 10^((BASE_POW+1) * 2) < 2^53, because of multiplyraw and divideraw
local BASE: number = 10 ^ BASE_POW

--> CONSTRUCTOR
local typeConstructors: {[string]: (any) -> ({number}, number)} = {
    ["string"] = function(str: string): ({number}, number)
        local isNegative: boolean = str:sub(1, 1) == "-"
        local strStart: number = isNegative and 2 or 1

        local arr: {number} = {}

        local i: number = #str
        repeat
            local chars: string = str:sub(max(i - BASE_POW + 1, strStart), i)
            local num: number? = tonumber(chars)
            if not num then break end

            i -= BASE_POW
            tinsert(arr, num)
        until i < strStart

        return arr, (isNegative and -1 or 1)
    end,

    ["number"] = function(num: number): ({number}, number)
        if num == 0 then
            return {0}, 0
        end

        local actualSign: number = sign(num)
        num = abs(num)

        local arr: {number} = {}
        while num > 0 do
            tinsert(arr, floor(num % BASE))
            num //= BASE
        end

        return arr, actualSign
    end,

    ["table"] = function(tbl: {any}): ({number}, number)
        local first: any = tbl[1]
        if not first then
            return {0}, 0
        end

        --> if is faster than assert, we use if here since its called a lot of times
        if typeof(first) ~= "number" then
            error("attempt to create AptInt from a table with non-int values")
        end

        return tbl, 1
    end,

    ["nil"] = function(_: nil): ({number}, number)
        return {0}, 0
    end,
}

-- =================
--> HELPER FUNCTIONS
-- =================


-- Split a digit into its low "d" limbs and its high limbs. For example, split_at(12345, 3) 
-- will extract the 3 final limbs, giving: high=12, low=345.
local function GetKaratsubaLower(int: AptInt, n: number): AptInt
    local limbs: {number} = int.limbs

    if #limbs <= n then
        return int:clone()
    end

    --> unsure if this is faster.
    if #limbs < SAFE_UNPACK_THRESHOLD then
        -- table.pack is faster here but FOR SOME REASON IT INCLUDES ["N"] IN THE ARRAY TOO?!
        return AptInt.new({unpack(limbs, 1, n)})
    end

    local result: {number} = tcreate(n, 0)

    for i: number = 1, n do
        result[i] = limbs[i]
    end

    return AptInt.new(result)
end

local function GetKaratsubaUpper(int: AptInt, n: number): AptInt
    -- {select(unpack(x))} is faster, but unpack only works up to a certain degree.
    local limbs: {number} = int.limbs
    local len: number = #limbs

    if len <= n then
        return AptInt.new()
    end

    if len < SAFE_UNPACK_THRESHOLD then
        -- select(n + 1, unpack(limbs)) is wayy faster than {unpack(limbs, n+1, #limbs)}
        --> for some reason we gotta set it to a variable first or else it breaks...
        local z: {number} = {select(n + 1, unpack(limbs))}
        return AptInt.new(z)
    end

    local result: {number} = {}

    for i: number = n + 1, len do
        tinsert(result, limbs[i])
    end

    return AptInt.new(result)
end


--> divsion stuff
local function CorrectRemainder(self: AptInt, divisor: AptInt, remainder: AptInt): AptInt
    local selfSign: number, divSign: number = self.signum, divisor.signum
    local sameSign: boolean = selfSign == divSign

    if sameSign then
        if divSign == -1 then
            return remainder:Negate(true)
        end

        return remainder
    elseif selfSign == 1 and divSign == -1 and remainder.signum ~= 0 then
        return divisor:AddRaw(remainder)
    elseif remainder.signum ~= 0 then
        return divisor:SubtractRaw(remainder)
    end

    return remainder
end

local function LowerThanAbsolute(x: AptInt, y: AptInt): boolean
    local xLimbs: {number}, yLimbs: {number} = x.limbs, y.limbs

    if #xLimbs < #yLimbs then
        return true
    end

    for i: number = #xLimbs, 1, -1 do
        local xDigit: number, yDigit: number = xLimbs[i], yLimbs[i] or 0
        if xDigit == yDigit then continue end

        return xDigit < yDigit
    end

    return false 
end

local function KaratsubaSqrt(n: AptInt): (AptInt, AptInt)
    --> NOTE: even tho we dont have burnikel-ziegler division, this is still faster by a ton
    
    -- https://gmplib.org/manual/Square-Root-Algorithm
    -- https://en.wikipedia.org/wiki/Integer_square_root#Karatsuba_square_root_algorithm
    if n.signum ~= 1 then
        return AptInt.new(), AptInt.new()
    end
    
    --> needed here too
    local limbLen: number = #n.limbs
    if limbLen < SQRT_KARATSUBA_THRESHOLD then
        local sq = n:sqrt()
        return sq, n:SubtractRaw(sq:MultiplyRaw(sq))
    end
    
    local AP_ONE: AptInt = AptInt.new({1})

    --> calculate a3, a2, a1 and a0
    local m2: number = limbLen // 2
    local m: number = m2 // 2
    
    local a3a2: AptInt, a1a0: AptInt = GetKaratsubaUpper(n, m2), GetKaratsubaLower(n, m2)
    local a1: AptInt, a0: AptInt = GetKaratsubaUpper(a1a0, m), GetKaratsubaLower(a1a0, m)
    
    local s1: AptInt, r1: AptInt = KaratsubaSqrt(a3a2)
    --local r1 = a3a2:SubtractRaw(s1:MultiplyRaw(s1))
    
    --> divrem (r1*b + a1, 2*s1)
    local q, u: AptInt = 
        r1:LeftShift(m):AddRaw(a1)
        :DivideRaw(s1:AddRaw(s1))

    --> get sqrt and correct
    local s: AptInt = s1:LeftShift(m):AddRaw(q) -- s=s1*b+q

    local r: AptInt = n:SubtractRaw(s:MultiplyRaw(s))
    
    if r.signum == -1 then
        r = r:AddRaw(s:AddRaw(s)):SubtractRaw(AP_ONE)
        s = s:SubtractRaw(AP_ONE)
    end

    return s, r
end

local function StripLeadingZeros(n: AptInt)
    local limbs: {number} = n.limbs
    while #limbs > 1 and limbs[#limbs] == 0 do
        limbs[#limbs] = nil
    end

    if #limbs == 1 and limbs[#limbs] == 0 then
        n.signum = 0
    end
end

-- ====================
--> METATABLE FUNCTIONS
-- ====================

-- Returns a new AptInt from the given argument.
function AptInt.new(num: ValidConstructor): AptInt
    local limbs: {number}, signum: number = typeConstructors[typeof(num)](num)

    return setmetatable({
        limbs = limbs,
        signum = signum
    }, AptInt) :: AptInt
end

--> this is ugly
local AP_NEG_ONE: AptInt = AptInt.new(-1)
local AP_ZERO: AptInt = AptInt.new()
local AP_ONE: AptInt = AptInt.new({1})
local AP_TWO: AptInt = AptInt.new({2})

-- Returns a copy of the AptInt.
function AptInt.clone(n: AptInt): AptInt
    return setmetatable({
        limbs = tclone(n.limbs),
        signum = n.signum
    }, AptInt) :: AptInt
end

-- Returns an AptInt with its sign flipped.
function AptInt:Negate(inPlace: boolean): AptInt
    if inPlace then
        self.signum = -self.signum
        return self
    end

    local result: AptInt = self:clone()
    result.signum = -result.signum

    return result
end

-- =====================
--> ARITHMETIC FUNCTIONS
-- =====================

-- Returns an AptInt equal to <code>a + b</code>.
function AptInt:AddRaw(term: AptInt): AptInt
    local selfSign: number, summandSign: number = self.signum, term.signum
    local sameSign: boolean = selfSign == summandSign

    --> CASES
    if summandSign == 0 then
        return self:clone()
    elseif selfSign == 0 then
        return term:clone()
    elseif selfSign == -1 and summandSign == 1 then -- -a+b = b+(-a) = b-a
        return (self:Negate(false):SubtractRaw(term)):Negate(true)
    elseif not sameSign then -- +a+(-b)=a-b
        return self:SubtractRaw(term:Negate(false))
    end

    local summandLimbs: {number}, selfLimbs: {number} = term.limbs, self.limbs
    if #selfLimbs < #summandLimbs then
        selfLimbs, summandLimbs =  summandLimbs, selfLimbs
    end

    local result: {number} = tcreate(#selfLimbs, 0)

    --> from right to left
    local carry: number = 0
    for i: number, selfDigit: number in selfLimbs do
        local sum: number = selfDigit + (summandLimbs[i] or 0) + carry

        if sum >= BASE then
            carry = 1
            sum = sum - BASE
        else
            carry = 0
        end

        result[i] = sum
    end

    if carry > 0 then
        result[#result + 1] = carry
    end

    local rAptInt: AptInt = AptInt.new(result)
    rAptInt.signum = (if sameSign and selfSign == -1 then -1 else 1)
    return rAptInt
end

-- Returns an AptInt, equal to <code>a - b</code>.
function AptInt:SubtractRaw(term: AptInt): AptInt
    local selfSign: number, termSign: number = self.signum, term.signum
    local sameSign: boolean = (selfSign == termSign)

    --> CASES
    if termSign == 0 then -- b - 0 = 0
        return self:clone()
    elseif selfSign == 0 then -- 0 - b = -b
        return term:Negate(false)
    elseif termSign == -1 and selfSign == 1 then -- a-(-b) = a+b
        return term:Negate(false):AddRaw(self)
    elseif sameSign and selfSign == -1 then -- (-a-(-b)) = b-a
        return (self:Negate(false):SubtractRaw(term:Negate(false))):Negate(true)
    end

    if self:LowerThanRaw(term) then
        return (term:SubtractRaw(self)):Negate(true)
    end

    local result: {number} = tcreate(#self.limbs, 0)

    local limbsSubrahend: {number} = term.limbs

    --> from right to left
    local borrow: number = 0
    for i: number, digit: number in self.limbs do
        local diff: number = digit - (limbsSubrahend[i] or 0) - borrow
        if diff < 0 then
            borrow = 1
            diff = diff + BASE
        else
            borrow = 0
        end

        result[i] = diff
    end

    --> StripLeadingZeros but on a table.
    while #result > 1 and result[#result] == 0 do
        result[#result] = nil
    end

    if #result == 1 and result[#result] == 0 then
        return AptInt.new()
    end

    return AptInt.new(result)
end

-- Returns an AptInt equal to <code>a * b</code>.
function AptInt:MultiplyRaw(factor: AptInt): AptInt
    local selfSign: number, factorSign: number = self.signum, factor.signum

    --> CASES
    if selfSign == 0 or factorSign == 0 then
        return AptInt.new()
    end

    local selflimbs: {number} = self.limbs
    local factorlimbs: {number} = factor.limbs
    local factorLen: number = #factorlimbs

    --> single limb multiplication
    if factorLen == 1 then
        local mul: number = factorlimbs[1]
        if mul == 1 then 
            local result = self:clone()
            result.signum = selfSign * factorSign

            return result
        end

        local result: AptInt = AptInt.new(tcreate(#selflimbs, 0))
        local rLimbs: {number} = result.limbs

        local carry: number = 0
        for i: number, limb: number in selflimbs do
            local p: number = limb * mul + carry

            carry = p // BASE
            rLimbs[i] = p % BASE
        end

        if carry > 0 then
            tinsert(rLimbs, carry)
        end

        result.signum = selfSign * factorSign
        return result
    end


    --> normal o(n^2) algorithm
    local len: number = #selflimbs + factorLen
    local result: AptInt = AptInt.new(tcreate(len, 0))
    local rLimbs: {number} = result.limbs

    local carry: number = 0
    for i: number = 1, len do
        local sum: number = carry
        for j: number, selfDigit: number in selflimbs do
            local index: number = i - j + 1
            if index < 1 or index > factorLen then continue end
            sum += selfDigit * factorlimbs[index]
        end

        rLimbs[i] = sum % BASE
        carry = sum // BASE
    end

    StripLeadingZeros(result)
    result.signum = selfSign * factorSign
    return result
end

-- Returns the quotient and remainder of a / b.
function AptInt:DivideRaw(divisor: AptInt): (AptInt, AptInt)
    --> we use knuths algorithm D for division

    local divSign: number, selfSign: number = divisor.signum, self.signum
    local sameSign: boolean = divSign == selfSign

    --> CASES (oh dear..)
    --> TODO: refactor
    if divSign == 0 or selfSign == 0 then -- dividing by 0
        return AptInt.new(), AptInt.new()
    elseif LowerThanAbsolute(self, divisor) then
        if selfSign == 1 and divSign == -1 then
            return AptInt.new({-1}), divisor:AddRaw(self)
        elseif selfSign == -1 and divSign == 1 then
            return AptInt.new({-1}), divisor:SubtractRaw(self:Negate(false))
        elseif sameSign and divSign == -1 then
            return AptInt.new(), self:clone()
        end

        return AptInt.new(), self:clone()
    elseif self:EqualsRaw(divisor) then -- dividing x/x = 1
        return AptInt.new({1}), AptInt.new()
    elseif self:Negate(false):EqualsRaw(divisor) or self:EqualsRaw(divisor:Negate(false)) then -- -x/x = -1
        return AptInt.new({-1}), AptInt.new()
    end

    local U: {number}, V: {number} = self.limbs, divisor.limbs

    --> single limb division
    if #V == 1 then
        local div: number = V[1]
        if div == 1 then
            local result = self:clone()
            result.signum = selfSign / divSign

            return result, AptInt.new()
        end

        local quotient: AptInt = AptInt.new(tcreate(#U - 1, 0))
        local qLimbs: {number} = quotient.limbs

        local carry: number = 0
        for i: number = #U, 1, -1 do
            local x: number = U[i] + carry * BASE
            local qDigit: number = x // div

            carry = x - qDigit * div
            qLimbs[i] = qDigit
        end

        quotient.signum = selfSign / divSign
        StripLeadingZeros(quotient)

        return quotient, AptInt.new({carry})
    end

    --> d0. defs
    local n: number = #V
    local m: number = #U - n

    --> d1. normalize (get Vn and Un)
    local D: number = BASE // (V[n] + 1)
    local Dint: AptInt = AptInt.new(D)

    local Vn: {number} = divisor:MultiplyRaw(Dint).limbs
    local Un: {number} = (self :: AptInt):MultiplyRaw(Dint).limbs

    local Vnn: number = Vn[n]
    local penultimateVn: number = Vn[n-1] or 0

    local quotient: AptInt = AptInt.new(tcreate(m + 1, 0))
    local qLimbs: {number} = quotient.limbs

    --> d2. init j
    local qhat: number, rhat: number = 0, 0
    for j: number = m, 0, -1 do
        --> d3. calculate qhat and rhat

        -- bigints for qhat and rhat are not needed, because un[x] is at most BASE
        -- currently, BASE * BASE + BASE < 2^53
        local nPlusJ: number = n + j + 1
        local top: number = (Un[nPlusJ] or 0) * BASE + (Un[nPlusJ - 1] or 0)
        qhat = top // Vnn
        rhat = top % Vnn

        if qhat == 0 then continue end

        --> first correction
        if qhat * penultimateVn > rhat * BASE + (Un[nPlusJ - 2] or 0) then
            qhat -= 1
            rhat += Vnn

            --> second correction
            if rhat < BASE and qhat * penultimateVn > rhat * BASE + (Un[nPlusJ - 2] or 0) then
                qhat -= 1
            end
        end

        --> d4. multiply & subtract
        local carry: number = 0
        local borrow: number = 0

        for i: number = 1, n do
            local p: number = qhat * Vn[i] + carry
            carry = p // BASE
            local pDigit: number = p % BASE

            local diff: number = Un[i + j] - pDigit - borrow
            if diff < 0 then
                borrow = 1
                Un[i + j] = diff + BASE
                continue
            end

            Un[i + j] = diff
            borrow = 0
        end

        local lastSub: number = (Un[nPlusJ] or 0) - carry - borrow
        Un[nPlusJ] = lastSub

        --> d6. add back
        if lastSub >= 0 then
            qLimbs[j + 1] = qhat
            continue
        end

        qLimbs[j + 1] = qhat - 1
        carry = 0

        for i: number = 1, n do
            local sum: number = Un[i + j] + Vn[i] + carry

            if sum >= BASE then
                carry = 1
                sum = sum - BASE
            else
                carry = 0
            end

            Un[i + j] = sum
        end

        Un[nPlusJ] = (Un[nPlusJ] + carry) % BASE
    end

    local remainder: AptInt = AptInt.new(Un)

    --> d8. unnormalize
    local carry: number = 0
    local rmLimbs: {number} = remainder.limbs
    for i: number = #rmLimbs, 1, -1 do
        local x: number = carry * BASE + rmLimbs[i]
        rmLimbs[i] = x // D :: number
        carry = x % D :: number
    end

    StripLeadingZeros(remainder)
    StripLeadingZeros(quotient)

    --> return
    remainder = CorrectRemainder(self, divisor, remainder)
    if not sameSign then
        quotient:Negate(true)
    end

    return quotient, remainder
end

-- Returns a new AptInt equal to sqrt(a).
function AptInt:sqrt(): AptInt
    --> we compute the isqrt via Newton-Heron iteration

    if self.signum ~= 1 then
        return AptInt.new()
    end
    
    if #(self :: AptInt).limbs >= SQRT_KARATSUBA_THRESHOLD then
        local s, _ = KaratsubaSqrt(self)
        return s
    end

    local prev2: AptInt = AP_NEG_ONE:clone()
    local prev1: AptInt = AP_ONE:LeftShift(#self.limbs//2+1)

    while true do
        local x1: AptInt = prev1:AddRaw(self:DivideRaw(prev1)):DivideRaw(AP_TWO)

        local equals: boolean = x1:EqualsRaw(prev1)
        if equals then
            return x1
        end

        if x1:EqualsRaw(prev2) and not equals then
            return (x1:LowerThanRaw(prev1) and x1 or prev1) -- min(x1, prev1)
        end

        prev2, prev1 = prev1, x1
    end
end

function AptInt:sqrt2(): AptInt
    --> we compute the isqrt via Newton-Heron iteration

    if self.signum ~= 1 then
        return AptInt.new()
    end

    local prev2: AptInt = AP_NEG_ONE:clone()
    local prev1: AptInt = AP_ONE:LeftShift(#self.limbs//2+1)

    while true do
        local x1: AptInt = prev1:AddRaw(self:DivideRaw(prev1)):DivideRaw(AP_TWO)

        local equals: boolean = x1:EqualsRaw(prev1)
        if equals then
            return x1
        end

        if x1:EqualsRaw(prev2) and not equals then
            return (x1:LowerThanRaw(prev1) and x1 or prev1) -- min(x1, prev1)
        end

        prev2, prev1 = prev1, x1
    end
end


-- ==================
--> QoL FUNCTIONS
-- ==================

-- Returns a new string representation of an AptInt
function AptInt:ToString(): string
    local str: string = (self.signum == -1 and "-" or "")
    if self.signum == 0 then
        return "0"
    end

    local limbs: {number} = self.limbs
    str ..=limbs[#limbs]

    local formatStr: string = `%0{BASE_POW}d`

    --> from left to right
    for i: number = #limbs - 1, 1, -1 do
        str ..= sformat(formatStr, limbs[i]) -- prepend dem zeros
    end

    return str
end

-- ==================
--> COMPARISON FUNCTIONS
-- ==================

-- Returns true if <code>self == num</code>.
function AptInt:EqualsRaw(num: AptInt): boolean
    local numLimbs: {number} = num.limbs
    if (#(self :: AptInt).limbs ~= #numLimbs) or (self.signum ~= num.signum) then
        return false
    end

    for i: number, selfDigit: number in self.limbs do
        if selfDigit ~= numLimbs[i] then
            return false
        end
    end

    return true
end

-- Returns true if <code>self < num</code>.
function AptInt:LowerThanRaw(num: AptInt): boolean
    local selfSign: number, numSign: number = self.signum, num.signum
    local selfLen: number, numLen: number = #self.limbs, #num.limbs

    --> check sign
    if selfSign < numSign then
        return true
    end

    local bothNegative: boolean = (selfSign == numSign and selfSign == -1)
    --> check their lengths. for negative, the lowest is with the most limbs, and for positive its the lowest
    if bothNegative and selfLen > numLen then
        return true
    elseif selfLen < numLen then -- both positive
        return true
    end

    local selflimbs: {number}, numlimbs: {number} = self.limbs, num.limbs
    for i: number = #selflimbs, 1, -1 do
        local sDigit: number, nDigit: number = selflimbs[i], numlimbs[i] or 0
        if sDigit == nDigit then continue end

        if bothNegative then
            return sDigit > nDigit
        end

        return sDigit < nDigit
    end

    return false 
end

-- Returns true if <code>self <= num</code>.
function AptInt:LowerOrEqualToRaw(num: AptInt): boolean
    local selfSign: number, numSign: number = self.signum, num.signum
    local selfLen: number, numLen: number = #self.limbs, #num.limbs

    --> check sign
    if selfSign <= numSign then
        return true
    end

    local bothNegative: boolean = (selfSign == numSign and selfSign == -1)
    --> check their lengths. for negative, the lowest is with the most limbs, and for positive its the lowest
    if bothNegative and selfLen >= numLen then
        return true
    elseif selfLen <= numLen then -- both positive
        return true
    end

    local selflimbs: {number}, numlimbs: {number} = self.limbs, num.limbs
    for i: number = #selflimbs, 1, -1 do
        local sDigit: number, nDigit: number = selflimbs[i], numlimbs[i] or 0
        if sDigit == nDigit then continue end

        if bothNegative then
            return sDigit >= nDigit
        end

        return sDigit <= nDigit
    end

    return true 
end

-- ==================
--> "BOOLEAN" FUNCTIONS
-- ==================

-- Returns a new AptInt equal to <code>n * 10<sup>(BASE * amount)</sup></code>
function AptInt:LeftShift(amount: number): AptInt
    local result: AptInt = self:clone()

    local tbl: {number} = result.limbs

    for i: number = 1, amount do
        tinsert(tbl, 1, 0)
    end

    return result
end

local apt1 = AptInt.new("143665816004337822710282600285310394341474369045835074863414468709543787931907367746179403311452034095731304066234112071267510472464260955153084575408147254672957261763907982395337943906645864229014250227057207826232751957053220218983971305018634078800548055251973907806245884614087189937340865371691338441989956445051526543084039211962387469415699218979531585795574920384684004258007709014706216763392717018544247025174258411677231986785008489302218244095")

-- IMPORTANT: to see if the sqrt was a success without any bugs, compare with :sqrt2() or with the bignumber calculator linked at the bottom.
-- IMPORTANT: karatsubasqrt only starts taking place at limbs >16*7 digits long (112 digits), because of the threshold

-- compare karatsubasqrt with base case newton-heron sqrt.
print(apt1:sqrt():EqualsRaw(apt1:sqrt2()))

-- expected: true
-- got: false

print(apt1:sqrt():ToString())
-- expected: 379032737378102767370356320425415662904513187772631008578870126471203845870697482014374611530431269030880793627229265919475483409207718357286202948008100864063587640630090308972232735749901964068667724412528434753635948938919935

-- got: 119860675788324263356489173083581170512663052701076259560030661381005856454449077924594509361740699212421677047855207038277359548991425890958622546606593188277707862040651192818989918644635073839432388646813517095663152474508

Above is the code needed to reproduce the bug. To see if Karatsubasqrt ran succesfully, compare :sqrt() with :sqrt2() (apt1:sqrt():EqualsRaw( apt1:sqrt2() ) . Karatsubasqrt also only kicks in at >121 digits.

Keep in mind that i am storing limbs using base 10^7 and LSD first. I've also tested all functions like DivideRaw and MultiplyRaw and they all work, its just my KaratsubaSqrt function thats wrong.

I'm using the base case function to compare my results (with the karatsuba condition removed), but I also double check using https://www.calculator.net/big-number-calculator.html.

Links to what I've been referencing for this are available at the top of the Karatsubasqrt code.

You can test code at https://luau.org/demo.

Answer 1

I have solved this. The problem lied in my splits, more specifically my m and m2 calculations, because they weren't even. The fixed code is:

--> split must be even
if limbLen % 2 ~= 0 then
    limbLen += 1
end
    
--> calculate a3, a2, a1 and a0
local m2: number = ceil(limbLen / 2)
if m2 % 2 ~= 0 then m2 += 1 end -- split must be even
local m: number = m2 / 2