Source code in the icon programming language

record num(val,err)

procedure main(a)
    x1 := num(100.0, 1.1)
    y1 := num(50.0,  1.2)
    x2 := num(200.0, 2.2)
    y2 := num(100.0, 2.3)
    d := pow(add(pow(sub(x1,x2),2),pow(sub(y1,y2),2)),0.5)
    write("d = [",d.val,", ",d.err,"]")
end

procedure add(a,b)
    return (numeric(a)+numeric(b)) |
           num(numeric(a)+b.val, b.err) |
           num(a.val+numeric(b), a.err) |
           num(a.val+b.val, (a.err^2 + b.err^2) ^ 0.5)
end

procedure sub(a,b)
    return (numeric(a)-numeric(b)) |
           num(numeric(a)-b.val, b.err) |
           num(a.val-numeric(b), a.err) |
           num(a.val-b.val, (a.err^2 + b.err^2) ^ 0.5)
end

procedure mul(a,b)
    return (numeric(a)*numeric(b)) |
           num(numeric(a)*b.val, abs(a*b.err)) |
           num(a.val*numeric(b), abs(b*a.err)) |
           num(f := a.val*b.val, ((f^2*((a.err/a.val)^2+(b.err/b.val)^2))^0.5))
end

procedure div(a,b)
    return (numeric(a)/numeric(b)) |
           num(numeric(a)/b.val, abs(a*b.err)) |
           num(a.val/numeric(b), abs(b*a.err)) |
           num(f := a.val/b.val, ((f^2*((a.err/a.val)^2+(b.err/b.val)^2))^0.5))
end

procedure pow(a,b)
    return (numeric(a)^numeric(b)) |
           num(f := a.val^numeric(b), abs(f*b*(a.err/a.val)))
end