Yash Raj Gupta / Jun 09 2019
Remix of

# Julia Season of Contributions Blog#1

My project titled 'Standard Compliant IntervalArithmetic Library' got selected in Julia Season of Contribution under the mentorship of David P. Sanders. My project is under JuliaIntervals to develop IntervalArithmetic.jl and make it compliant with the IEEE1788-2015, a paper which provides standards for interval computations. Basically my project comprised of including the missing functions and improving the existing functions upto the standards.

## Community Bonding

I started with reading the standard for the implementation of the Reduction and String Parsing functions. In the meantime, I also started to solve some issues, namely `cospi/sinpi stack overflow` and `broken rational power functions`. I learned how to write regular expressions for the string parsing functions.

## Starting with the project

During the end of the community boding period I started with the implementation of the functions planned in my proposal. I had reduction functions and the string parsing functions in my first two weeks. Lets see how did they go in the implementation section.

## Implementation

### Sinpi and Cospi functions

`sinPi` and `cosPi` are defined as sin(πx) and cos(πx) respectively with domain as R and range equal to [-1, 1]. These functions should return a range containing the sine and cosine of values in the interval πx. I used the implementation of sin and cos functions to implement `sinPi` and `cosPi` functions in IntervalArithmetic.jl.

```import Pkg
using IntervalArithmetic
function sinpi(a::Interval{T}) where T
isempty(a) && return a
w = a * pi_interval(T)
return sin(w)
end
sinpi(0.5 .. 0.75)```
[0.707106, 1]

### Fixing Rational Power

The function `^(a :: Interval, x :: Rational)` computes the range of values in the interval raised to the power x. This issue was not straight as it was ambiguous on some platform. First I started with using the power functions by type casting `x `into BigFloat. The added tests worked well but it failed for some existing tests. Then I moved onto compute the qth root raised to the power p where `p` is numerator and `q` is the denominator of `x`. I calculated the lower and upper bounds for qth root but I could not get the correct rounding as per the standard. Then I used `ccall` function to wrap `nthroot` function in GNU MPFR library to ensure correct rounded results.

```function ^(a::Interval{T}, x::Rational) where T
if x > 0
low = a.lo ^ BigFloat(x)
high = a.hi ^ BigFloat(x)
if isinteger(high) && isinteger(low)
return Interval(low, high)
end
return ^(a, BigFloat(x))
end
if x < 0
return inv(^(a, -x))
end
end```
^ (generic function with 1 method)
```function ^(a::Interval{T}, x::Rational) where T
p = x.num
q = x.den
if x > 0
low = a.lo ^ BigFloat(1//q)
high = a.hi ^ BigFloat(1//q)
isinteger(high) && isinteger(low) && return ^(Interval(low, high) , p)
b = @interval(low, high)
return b ^ p
end
```
0.6s
```import Base.MPFR: MPFRRoundingMode
import Base.MPFR: MPFRRoundUp, MPFRRoundDown, MPFRRoundNearest, MPFRRoundToZero, MPFRRoundFromZero
function ^(a::Interval{T}, x::Rational) where T
domain = Interval{T}(0, Inf)

p = x.num
q = x.den

isempty(a) && return emptyinterval(a)
x == 0 && return one(a)
if a == zero(a)
x > zero(x) && return zero(a)
return emptyinterval(a)
end
x == (1//2) && return sqrt(a)

if x >= 0
if a.lo ≥ 0
isinteger(x) && return a ^ Int64(x)
a = @biginterval(a)
ui = convert(Culong, q)
low = BigFloat()
high = BigFloat()
ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown)
ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp)
b = @interval(low , high)
return b^p
end
if a.lo < 0 && a.hi ≥ 0
isinteger(x) && return a ^ Int64(x)
a = a ∩ Interval{T}(0, Inf)
a = @biginterval(a)
ui = convert(Culong, q)
low = BigFloat()
high = BigFloat()
ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown)
ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp)
b = @interval(low , high)
return b^p
end

if a.hi < 0
isinteger(x) && return a ^ Int64(x)
return emptyinterval(a)
end
end

if x < 0
return inv(a^(-x))
end
end```
^ (generic function with 1 method)

### Reduction functions

We have four operations on vectors, namely sum, dot, sumAbs and sumSquare. This implementation is currently under review. There are issues related to the precision of BigFloat and generating different tests for different rounding modes. The standard states that these functions should return the correctly rounded results hence I used `BigFloat` and then rounded it to return `Float64` value.

```function vector_sum(v :: Vector{T} , r :: RoundingMode) where T
sum1 = sum(BigFloat.(v))
return (sum1, r)
end

function vector_dot(v :: Vector{T}, u :: Vector{T}, r :: RoundingMode) where T
sum1 = sum(BigFloat.(v) .* BigFloat.(u))
return (sum1, r)
end

function vector_sumAbs(v :: Vector{T}, r :: RoundingMode) where T
sum1 = sum(BigFloat.(abs.(v)))
return (sum1, r)
end

function vector_sumSquare(v :: Vector{T}, r :: RoundingMode) where T
return vector_dot(v, v, r)
end```
(15.0, RoundingMode{:Up}())

### String Parsing

I had to implement mostly the uncertain form(like 2.4?1) conversion. Moreover it required recognition of regular expressions for different types of strings.

0.6s
```function parse(::Type{Interval{T}}, s::AbstractString) where T
m = match(r"(\-?\d*\.?\d*)\?\?([ud]?)", s)   # match with string of form "34??"

if m!= nothing
if m.captures[2] == ""    # strings of the form "10??"
return entireinterval(T)
end
if m.captures[2] == "u"   # strings of the form "10??u"
lo = parse(Float64, m.captures[1])
hi = Inf
interval = eval(make_interval(T, lo, [hi]))
return interval
end
if m.captures[2] == "d"   # strings of the form "10??d"
lo = -Inf
hi = parse(Float64, m.captures[1])
interval = eval(make_interval(T, lo, [hi]))
return interval
end
end

m = match(r"(\-?\d*)\.?(\d*)\?(\d*)([ud]?)(e-?\d*)?", s) # match with strings like "3.56?1u" or "3.56?1e2"

if m != nothing  # matched
if m.captures[3] != "" && m.captures[4] == "" && m.captures[5] == nothing # string of form "3.4?1" or "10?2"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = n - x
hi = n + x
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] == nothing  # strings of the form "3.46?"
d = length(m.captures[2])
x = parse(Float64, "0.5" * "e-\$d")
n = parse(Float64, m.captures[1]*"." * m.captures[2])
lo = n - x
hi = n + x
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] != nothing  # strings of the form "3.46?e2"
d = length(m.captures[2])
x = parse(Float64, "0.5" * "e-\$d")
n = parse(Float64, m.captures[1]*"." * m.captures[2])
lo = parse(Float64, string(n-x) * m.captures[5])
hi = parse(Float64, string(n+x) * m.captures[5])
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[3] == "" && m.captures[4] == "u" && m.captures[5] == nothing # strings of the form "3.4?u"
d = length(m.captures[2])
x = parse(Float64, "0.5" * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = n
hi = n+x
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[3] == "" && m.captures[4] == "d" && m.captures[5] == nothing # strings of the form "3.4?d"
d = length(m.captures[2])
x = parse(Float64, "0.5" * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = n - x
hi = n
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[3] != ""
if m.captures[4] == "u" && m.captures[5] != nothing    #strings of the form "3.46?1u"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = parse(Float64, string(n) * m.captures[5])
hi = parse(Float64, string(n+x) * m.captures[5])
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[4] == "u" && m.captures[5] == nothing    #strings of the form "3.46?1u"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = n
hi = n+x
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[4] == "d" && m.captures[5] != nothing  #strings of the form "3.46?1d"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = parse(Float64, string(n-x) * m.captures[5])
hi = parse(Float64, string(n) * m.captures[5])
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[4] == "d" && m.captures[5] == nothing  #strings of the form "3.46?1d"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = n-x
hi = n
interval = eval(make_interval(T, lo, [hi]))
return interval
end

if m.captures[4] == "" && m.captures[5] != nothing    # strings of the form "3.56?1e2"
d = length(m.captures[2])
x = parse(Float64, m.captures[3] * "e-\$d")
n = parse(Float64, m.captures[1]*"."*m.captures[2])
lo = parse(Float64, string(n-x)*m.captures[5])
hi = parse(Float64, string(n+x)*m.captures[5])
interval = eval(make_interval(T, lo, [hi]))
return interval
end
end

m = match(r"(-?\d*\.?\d*)", s)  # match strings of form "1" and "2.4"

if m != nothing
lo = parse(Float64, m.captures[1])
hi = lo
end

if m == nothing
throw(ArgumentError("Unable to process string \$s as interval"))
end

interval = eval(make_interval(T, lo, [hi]))
return interval
end
end```
parse (generic function with 1 method)

## What next?

I am following my project timeline and will try to do so during whole project. I have IntervalToText conversion next on my list and some other functions. I will also try to fix some issues apart from these. I must say its been great experience with the Julia community till now and I am looking forward to the completion of my project.