Arithmetics: Implement Base.sqrt (#102)

Author: barucden

docs/src/operations.md

@@ -89,6 +89,19 @@ julia> dec"0" / dec"0"
 ERROR: UndefinedDivisionError()
 ```
 
+### Square root
+
+![Affected by context](https://img.shields.io/badge/ctxt-affected-blue)
+
+Square root is implemented via the `Base.sqrt` function.
+```jldoctest
+julia> sqrt(dec"9")
+3
+
+julia> sqrt(dec"2")
+1.414213562373095048801688724
+```
+
 ### Absolute value
 
 ![Affected by context](https://img.shields.io/badge/ctxt-affected-blue)

scripts/dectest.jl

@@ -47,6 +47,7 @@ function translate(io, dectest_path)
 
             test = parse_test(line)
             any(isspecial, test.operands) && continue
+            isspecial(test.result) && continue
 
             dectest = decimal_test(test, directives)
             println(io, dectest)
@@ -60,7 +61,7 @@ function isspecial(value)
 end
 
 function parse_precision(line)
-    m = match(r"^precision:\s*(\d+)$", line)
+    m = match(r"^precision:\s*(\d+).*$", line)
     isnothing(m) && throw(ArgumentError(line))
     return parse(Int, m[1])
 end
@@ -187,6 +188,8 @@ function decimal_operation(operation, operands)
         return decimal_reduce(operands...)
     elseif operation == "subtract"
         return decimal_subtract(operands...)
+    elseif operation == "squareroot"
+        return decimal_sqrt(operands...)
     else
         throw(ArgumentError(operation))
     end
@@ -204,4 +207,5 @@ decimal_multiply(x, y) = :($(dec(x)) * $(dec(y)))
 decimal_plus(x) = :(+($(dec(x))))
 decimal_reduce(x) = :(normalize($(dec(x))))
 decimal_subtract(x, y) = :($(dec(x)) - $(dec(y)))
+decimal_sqrt(x) = :(sqrt($(dec(x))))
 

src/Decimals.jl

@@ -5,7 +5,6 @@
 module Decimals
 
 export Decimal,
-       number,
        normalize,
        @dec_str,
        DivisionByZeroError,
@@ -26,7 +25,9 @@ include("bigint.jl")
 include("context.jl")
 include("conversion.jl")
 include("decimal.jl")
-include("arithmetic.jl")
+include("arithmetic/elementary.jl")
+include("arithmetic/exceptions.jl")
+include("arithmetic/sqrt.jl")
 include("equals.jl")
 include("round.jl")
 include("hash.jl")

src/arithmetic.jl src/arithmetic/elementary.jlsrc/arithmetic.jl renamed to src/arithmetic/elementary.jl

@@ -1,23 +1,7 @@
 Base.promote_rule(::Type{Decimal}, ::Type{<:Real}) = Decimal
-
-# override definitions in Base
 Base.promote_rule(::Type{BigFloat}, ::Type{Decimal}) = Decimal
 Base.promote_rule(::Type{BigInt}, ::Type{Decimal}) = Decimal
 
-"""
-    DivisionByZeroError
-
-Division was attempted with a denominator value of 0.
-"""
-struct DivisionByZeroError <: Exception end
-
-"""
-    UndefinedDivisionError
-
-Division was attempted with both numerator and denominator value of 0.
-"""
-struct UndefinedDivisionError <: Exception end
-
 Base.:(+)(x::Decimal) = fix(x)
 Base.:(-)(x::Decimal) = fix(Decimal(!x.s, x.c, x.q))
 
@@ -220,4 +204,3 @@ end
 
 Base.abs(x::Decimal) = fix(Decimal(false, x.c, x.q))
 
-# TODO exponentiation

src/arithmetic/exceptions.jl

@@ -0,0 +1,14 @@
+"""
+    DivisionByZeroError
+
+Division was attempted with a denominator value of 0.
+"""
+struct DivisionByZeroError <: Exception end
+
+"""
+    UndefinedDivisionError
+
+Division was attempted with both numerator and denominator value of 0.
+"""
+struct UndefinedDivisionError <: Exception end
+

src/arithmetic/sqrt.jl

@@ -0,0 +1,57 @@
+function Base.sqrt(x::Decimal)
+    if iszero(x)
+        return x
+    end
+
+    if signbit(x)
+        throw(DomainError(x, "Square root of decimals is defined only for non-negative arguments"))
+    end
+
+    if isone(x)
+        return x
+    end
+
+    (; c, q) = x
+
+    # We are computing
+    #
+    #   sqrt(c * 10^q) = sqrt(c) * sqrt(10^q)
+    #                  =      d  * sqrt(10^q)
+    #
+    # where `d` is constrained to have `precision` digits:
+    #
+    #         10^(p - 1) ≤       d ≤ 10^p
+    #         10^(p - 1) ≤ sqrt(c) ≤ 10^p
+    #   10^(2 * (p - 1)) ≤       c ≤ 10^(2p)
+    #
+    # Consequently, we need `c` to have at least `2 * (precision - 1)` digits.
+    # However, to figure out correct rounding, we compute in an increased
+    # precision, `precision + 1`.
+    n = 1 + 2 * precision(Decimal) - ndigits(c)
+    if n > 0
+        c = c * BigTen^n
+        q = q - n
+    end
+
+    # To make sure, that we can compute `sqrt(10^q)`, we make sure that `q` is
+    # even, `q = 2r` so that `sqrt(10^q) = sqrt(10^2r) = 10^r`.
+    if isodd(q)
+        c = 10 * c
+        q -= 1
+    end
+
+    d = isqrt(c)
+    r = q ÷ 2
+
+    # If `d % 5 == 0`, we add one to `d` so that `fix` takes care of the
+    # rounding properly
+    if isdivisible(d, 5) && d^2 != c
+        d += 1
+    end
+
+    d, m = cancelfactor(d, Val(10))
+    r += m
+
+    return fix(Decimal(0, d, r))
+end
+

test/arithmetic/test_sqrt.jl

@@ -0,0 +1,9 @@
+using Decimals
+using Test
+
+@testset "sqrt" begin
+    @test iszero(sqrt(zero(Decimal)))
+    @test isone(sqrt(one(Decimal)))
+    @test_throws DomainError sqrt(Decimal(1, 1, 0))
+end
+