Numeric error propagation

From Rosetta Code
Jump to: navigation, search
Task
Numeric error propagation
You are encouraged to solve this task according to the task description, using any language you may know.

If f, a, and b are values with uncertainties σf, σa, and σb. and c is a constant; then if f is derived from a, b, and c in the following ways, then σf can be calculated as follows:

Addition/Subtraction
  • If f = a ± c, or f = c ± a then σf = σa
  • If f = a ± b then σf2 = σa2 + σb2
Multiplication/Division
  • If f = ca or f = ac then σf = |cσa|
  • If f = ab or f = a / b then σf2 = f2( (σa / a)2 + (σb / b)2)
Exponentiation
  • If f = ac then σf = |fc(σa / a)|
Caution:
This implementation of error propagation does not address issues of dependent and independent values. It is assumed that a and b are independent and so the formula for multiplication should not be applied to a*a for example. See the talk page for some of the implications of this issue.
Task details
  1. Add an uncertain number type to your language that can support addition, subtraction, multiplication, division, and exponentiation between numbers with an associated error term together with 'normal' floating point numbers without an associated error term.
    Implement enough functionality to perform the following calculations.
  2. Given coordinates and their errors:
    x1 = 100 ± 1.1
    y1 = 50 ± 1.2
    x2 = 200 ± 2.2
    y2 = 100 ± 2.3
    if point p1 is located at (x1, y1) and p2 is at (x2, y2); calculate the distance between the two points using the classic pythagorean formula:
    d = √((x1 - x2)2 + (y1 - y2)2)
  3. Print and display both d and its error.
References
Cf.

Contents

[edit] Ada

Specification of a generic type Approximation.Number, providing all the operations required to solve the task ... and some more operations, for completeness.

generic
type Real is digits <>;
with function Sqrt(X: Real) return Real;
with function "**"(X: Real; Y: Real) return Real;
package Approximation is
 
type Number is private;
 
-- create an approximation
function Approx(Value: Real; Sigma: Real) return Number;
 
-- unary operations and conversion Real to Number
function "+"(X: Real) return Number;
function "-"(X: Real) return Number;
function "+"(X: Number) return Number;
function "-"(X: Number) return Number;
 
-- addition / subtraction
function "+"(X: Number; Y: Number) return Number;
function "-"(X: Number; Y: Number) return Number;
 
-- multiplication / division
function "*"(X: Number; Y: Number) return Number;
function "/"(X: Number; Y: Number) return Number;
 
-- exponentiation
function "**"(X: Number; Y: Positive) return Number;
function "**"(X: Number; Y: Real) return Number;
 
-- Output to Standard IO (wrapper for Ada.Text_IO and Ada.Text_IO.Float_IO)
procedure Put_Line(Message: String;
Item: Number;
Value_Fore: Natural := 7;
Sigma_Fore: Natural := 4;
Aft: Natural := 2;
Exp: Natural := 0);
procedure Put(Item: Number;
Value_Fore: Natural := 7;
Sigma_Fore: Natural := 3;
Aft: Natural := 2;
Exp: Natural := 0);
 
private
type Number is record
Value: Real;
Sigma: Real;
end record;
end Approximation;

The implementation:

with Ada.Text_IO;
 
package body Approximation is
 
package RIO is new Ada.Text_IO.Float_IO(Real);
 
-- create an approximation
 
function Approx(Value: Real; Sigma: Real) return Number is
begin
return (Value, Sigma);
end Approx;
 
-- unary operations and conversion Real to Number
 
function "+"(X: Real) return Number is
begin
return Approx(X, 0.0);
end "+";
 
function "-"(X: Real) return Number is
begin
return Approx(-X, 0.0);
end "-";
 
function "+"(X: Number) return Number is
begin
return X;
end "+";
 
function "-"(X: Number) return Number is
begin
return Approx(-X.Value, X.Sigma);
end "-";
 
-- addition / subtraction
 
function "+"(X: Number; Y: Number) return Number is
Z: Number;
begin
Z.Value := X.Value + Y.Value;
Z.Sigma := Sqrt(X.Sigma*X.Sigma + Y.Sigma*Y.Sigma);
return Z;
end "+";
 
function "-"(X: Number; Y: Number) return Number is
begin
return X + (-Y);
end "-";
 
-- multiplication / division
 
function "*"(X: Number; Y: Number) return Number is
Z: Number;
begin
Z.Value := X.Value * Y.Value;
Z.Sigma := Z.Value * Sqrt((X.Sigma/X.Value)**2 + (Y.Sigma/Y.Value)**2);
return Z;
end "*";
 
function "/"(X: Number; Y: Number) return Number is
Z: Number;
begin
Z.Value := X.Value / Y.Value;
Z.Sigma := Z.Value * Sqrt((X.Sigma/X.Value)**2 + (Y.Sigma/Y.Value)**2);
return Z;
end "/";
 
-- exponentiation
 
function "**"(X: Number; Y: Positive) return Number is
Z: Number;
begin
Z.Value := X.Value ** Y ;
Z.Sigma := Z.Value * Real(Y) * (X.Sigma/X.Value);
if Z.Sigma < 0.0 then
Z.Sigma := - Z.Sigma;
end if;
return Z;
end "**";
 
function "**"(X: Number; Y: Real) return Number is
Z: Number;
begin
Z.Value := X.Value ** Y ;
Z.Sigma := Z.Value * Y * (X.Sigma/X.Value);
if Z.Sigma < 0.0 then
Z.Sigma := - Z.Sigma;
end if;
return Z;
end "**";
 
-- Output to Standard IO (wrapper for Ada.Text_IO.Float_IO)
 
procedure Put_Line(Message: String;
Item: Number;
Value_Fore: Natural := 7;
Sigma_Fore: Natural := 4;
Aft: Natural := 2;
Exp: Natural := 0) is
begin
Ada.Text_IO.Put(Message);
Put(Item, Value_Fore, Sigma_Fore, Aft, Exp);
Ada.Text_IO.New_Line;
end Put_Line;
 
procedure Put(Item: Number;
Value_Fore: Natural := 7;
Sigma_Fore: Natural := 3;
Aft: Natural := 2;
Exp: Natural := 0) is
begin
RIO.Put(Item.Value, Value_Fore, Aft, Exp);
Ada.Text_IO.Put(" (+-");
RIO.Put(Item.Sigma, Sigma_Fore, Aft, Exp);
Ada.Text_IO.Put(")");
end Put;
 
end Approximation;

Instantiating the package with Float operations, to compute the distance:

with Approximation, Ada.Numerics.Elementary_Functions;
 
procedure Test_Approximations is
package A is new Approximation(Float,
Ada.Numerics.Elementary_Functions.Sqrt,
Ada.Numerics.Elementary_Functions."**");
use type A.Number;
X1: A.Number := A.Approx(100.0, 1.1);
Y1: A.Number := A.Approx( 50.0, 1.2);
X2: A.Number := A.Approx(200.0, 2.2);
Y2: A.Number := A.Approx(100.0, 2.3);
 
begin
A.Put_Line("Distance:",
((X1-X2)**2 + (Y1 - Y2)**2)**0.5,
Sigma_Fore => 1);
end Test_Approximations;

Output:

Distance:    111.80 (+-2.49)

[edit] C

Rewrote code to make it more compact and added a nice formatting function for imprecise values so that they are printed out in a technically correct way i.e. with the symbol '±' . Output pasted after code.

 
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
 
typedef struct{
double value;
double delta;
}imprecise;
 
#define SQR(x) ((x) * (x))
imprecise imprecise_add(imprecise a, imprecise b)
{
imprecise ret;
ret.value = a.value + b.value;
ret.delta = sqrt(SQR(a.delta) + SQR(b.delta));
return ret;
}
 
imprecise imprecise_mul(imprecise a, imprecise b)
{
imprecise ret;
ret.value = a.value * b.value;
ret.delta = sqrt(SQR(a.value * b.delta) + SQR(b.value * a.delta));
return ret;
}
 
imprecise imprecise_div(imprecise a, imprecise b)
{
imprecise ret;
ret.value = a.value / b.value;
ret.delta = sqrt(SQR(a.value * b.delta) + SQR(b.value * a.delta)) / SQR(b.value);
return ret;
}
 
imprecise imprecise_pow(imprecise a, double c)
{
imprecise ret;
ret.value = pow(a.value, c);
ret.delta = fabs(ret.value * c * a.delta / a.value);
return ret;
}
 
char* printImprecise(imprecise val)
{
char principal[30],error[30],*string,sign[2];
sign[0] = 241; /* ASCII code for ±, technical notation for denoting errors */
sign[1] = 00;
 
sprintf(principal,"%f",val.value);
sprintf(error,"%f",val.delta);
 
string = (char*)malloc((strlen(principal)+1+strlen(error)+1)*sizeof(char));
 
strcpy(string,principal);
strcat(string,sign);
strcat(string,error);
 
return string;
}
 
int main(void) {
imprecise x1 = {100, 1.1};
imprecise y1 = {50, 1.2};
imprecise x2 = {-200, 2.2};
imprecise y2 = {-100, 2.3};
imprecise d;
 
d = imprecise_pow(imprecise_add(imprecise_pow(imprecise_add(x1, x2), 2),imprecise_pow(imprecise_add(y1, y2), 2)), 0.5);
printf("Distance, d, between the following points :");
printf("\n( x1, y1) = ( %s, %s)",printImprecise(x1),printImprecise(y1));
printf("\n( x2, y2) = ( %s, %s)",printImprecise(x2),printImprecise(y2));
printf("\nis d = %s", printImprecise(d));
return 0;
}
 
Distance, d, between the following points :
( x1, y1) = ( 100.000000±1.100000, 50.000000±1.200000)
( x2, y2) = ( -200.000000±2.200000, -100.000000±2.300000)
is d = 111.803399±2.487167

[edit] D

import std.stdio, std.math, std.string, std.typecons, std.traits;
 
const struct Imprecise {
private const double value, delta;
 
this(in double v, in double d) pure nothrow {
this.value = v;
this.delta = abs(d);
}
 
enum IsImprecise(T) = is(Unqual!T == Unqual!(typeof(this)));
 
I reciprocal() const pure nothrow {
return I(1.0 / value, delta / (value ^^ 2));
}
 
string toString() const {
return format("I(value=%g, delta=%g)", value, delta);
}
 
I opUnary(string op:"-")() const pure nothrow {
return I(-this.value, this.delta);
}
 
I opBinary(string op:"+", T)(in T other) const pure nothrow
if (isNumeric!T || IsImprecise!T) {
static if (IsImprecise!T)
return I(this.value + other.value,
(this.delta ^^ 2 + other.delta ^^ 2) ^^ 0.5);
else
return I(this.value + other, this.delta);
}
 
I opBinaryRight(string op:"+", T)(in T other) const pure nothrow
if (isNumeric!T) {
return I(this.value + other, this.delta);
}
 
I opBinary(string op:"-", T)(in T other) const pure nothrow
if (isNumeric!T || IsImprecise!T) {
return this + (-other);
}
 
I opBinaryRight(string op:"-", T)(in T other) const pure nothrow
if (isNumeric!T) {
return this - other;
}
 
I opBinary(string op:"*", T)(in T other) const pure nothrow
if (isNumeric!T || IsImprecise!T) {
static if (IsImprecise!T) {
auto f = this.value * other.value;
return I(f, f * ((delta / value) ^^ 2 +
(other.delta / other.value) ^^ 2) ^^ 0.5);
} else
return I(this.value * other, this.delta * other);
}
 
I opBinaryRight(string op:"*", T)(in T other) const pure nothrow
if (isNumeric!T) {
return this * other;
}
 
I opBinary(string op:"/", T)(in T other) const pure nothrow
if (isNumeric!T || IsImprecise!T) {
static if (IsImprecise!T)
return this * other.reciprocal();
else
return I(this.value / other, this.delta / other);
}
 
I opBinaryRight(string op:"/", T)(in T other) const pure nothrow
if (isNumeric!T) {
return this / other;
}
 
I opBinary(string op:"^^", T)(in T other) const pure nothrow
if (isNumeric!T) {
auto f = this.value ^^ other;
return I(f, f * other * (this.delta / this.value));
}
}
 
alias I = Imprecise;
 
auto distance(T1, T2)(in T1 p1, in T2 p2) pure nothrow {
return ((p1[0] - p2[0]) ^^ 2 + (p1[1] - p2[1]) ^^ 2) ^^ 0.5;
}
 
void main() {
immutable x1 = I(100, 1.1);
immutable x2 = I(200, 2.2);
immutable y1 = I( 50, 1.2);
immutable y2 = I(100, 2.3);
 
immutable p1 = tuple(x1, y1);
immutable p2 = tuple(x2, y2);
writefln("Point p1: (%s, %s)", p1[0], p1[1]);
writefln("Point p2: (%s, %s)", p2[0], p2[1]);
writeln("Distance(p1, p2): ", distance(p1, p2));
}
Output:
Point p1: (I(value=100, delta=1.1), I(value=50, delta=1.2))
Point p2: (I(value=200, delta=2.2), I(value=100, delta=2.3))
Distance(p1, p2): I(value=111.803, delta=2.48717)

[edit] Go

Variance from task requirements is that the following does not "extend the language." It simply defines a type with associated functions and methods as required to solve the remainder of the task.

package main
 
import (
"fmt"
"math"
)
 
// "uncertain number type"
// a little optimization is to represent the error term with its square.
// this saves some taking of square roots in various places.
type unc struct {
n float64 // the number
s float64 // *square* of one sigma error term
}
 
// constructor, nice to have so it can handle squaring of error term
func newUnc(n, s float64) *unc {
return &unc{n, s * s}
}
 
// error term accessor method, nice to have so it can handle recovering
// (non-squared) error term from internal (squared) representation
func (z *unc) errorTerm() float64 {
return math.Sqrt(z.s)
}
 
// Arithmetic methods are modeled on the Go big number package.
// The basic scheme is to pass all operands as method arguments, compute
// the result into the method receiver, and then return the receiver as
// the result of the method. This has an advantage of letting the programer
// determine allocation and use of temporary objects, reducing garbage;
// and has the convenience and efficiency of allowing operations to be chained.
 
// addition/subtraction
func (z *unc) addC(a *unc, c float64) *unc {
*z = *a
z.n += c
return z
}
 
func (z *unc) subC(a *unc, c float64) *unc {
*z = *a
z.n -= c
return z
}
 
func (z *unc) addU(a, b *unc) *unc {
z.n = a.n + b.n
z.s = a.s + b.s
return z
}
func (z *unc) subU(a, b *unc) *unc {
z.n = a.n - b.n
z.s = a.s + b.s
return z
}
 
// multiplication/division
func (z *unc) mulC(a *unc, c float64) *unc {
z.n = a.n * c
z.s = a.s * c * c
return z
}
 
func (z *unc) divC(a *unc, c float64) *unc {
z.n = a.n / c
z.s = a.s / (c * c)
return z
}
 
func (z *unc) mulU(a, b *unc) *unc {
prod := a.n * b.n
z.n, z.s = prod, prod*prod*(a.s/(a.n*a.n)+b.s/(b.n*b.n))
return z
}
 
func (z *unc) divU(a, b *unc) *unc {
quot := a.n / b.n
z.n, z.s = quot, quot*quot*(a.s/(a.n*a.n)+b.s/(b.n*b.n))
return z
}
 
// exponentiation
func (z *unc) expC(a *unc, c float64) *unc {
f := math.Pow(a.n, c)
g := f * c / a.n
z.n = f
z.s = a.s * g * g
return z
}
 
func main() {
x1 := newUnc(100, 1.1)
x2 := newUnc(200, 2.2)
y1 := newUnc(50, 1.2)
y2 := newUnc(100, 2.3)
var d, d2 unc
d.expC(d.addU(d.expC(d.subU(x1, x2), 2), d2.expC(d2.subU(y1, y2), 2)), .5)
fmt.Println("d: ", d.n)
fmt.Println("error:", d.errorTerm())
}

Output:

d:     111.80339887498948
error: 2.487167063146342

[edit] Haskell

data Error a = Error {value :: a, uncertainty :: a} deriving (Eq, Show)
 
instance (Floating a) => Num (Error a) where
Error a ua + Error b ub = Error (a + b) (sqrt (ua ^ 2 + ub ^ 2))
negate (Error a ua) = Error (negate a) ua
Error a ua * Error b ub = Error (a * b) (abs (a * b * sqrt ((ua / a) ^ 2 + (ub / b) ^ 2))) -- I've factored out the f^2 from the square root
fromInteger a = Error (fromInteger a) 0
 
instance (Floating a) => Fractional (Error a) where
fromRational a = Error (fromRational a) 0
Error a ua / Error b ub = Error (a / b) (abs (a / b * sqrt ((ua / a) ^ 2 + (ub / b) ^ 2))) -- I've factored out the f^2 from the square root
 
instance (Floating a) => Floating (Error a) where
Error a ua ** Error c 0 = Error (a ** c) (abs (ua * c * a**c / a))
 
main = print (sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2)) where -- using (^) for exponentiation would calculate a*a, which the problem specifically said was calculated wrong
x1 = Error 100 1.1
y1 = Error 50 1.2
x2 = Error 200 2.2
y2 = Error 100 2.3
 
Output:
Error {value = 111.80339887498948, uncertainty = 2.4871670631463423}

[edit] Icon and Unicon

The following solution works in both languages.

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

The output is:

->nep
d = [111.8033988749895, 2.487167063146342]
->

[edit] J

J's built in operators cannot be overloaded to deal with user defined types. So we will have to create new operators. Here's one approach, which is sufficient for this example:

First, we will need some utilities. num will extract the number part of a number. unc will extract the uncertainty part of a number, and will also be used to associate uncertainty with a number. dist will compute the distance between two numbers (which is needed for multiplicative uncertainty).

num=: {."1
unc=: {:@}."1 : ,.
dist=: +/&.:*:

Jumping into the example values, for illustration purposes:

x1=: 100 unc 1.1
y1=: 50 unc 1.2
 
x2=: 200 unc 2.2
y2=: 100 unc 2.3

Above, we see unc being used to associate a number with its uncertainty. Here's how to take them apart again:

   num x1
100
unc x1
1.1

Note that these operations "do the right thing" for normal numbers:

   num 100
100
unc 100
0

And, a quick illustration of the distance function:

   3 dist 4
5

Next, we need to define our arithmetic operations:

add=: +&num unc dist&unc
sub=: -&num unc dist&unc
mul=: *&num unc |@(*&num * dist&(unc%num))
div=: %&num unc |@(%&num * dist&(unc%num))
exp=: ^&num unc |@(^&num * dist&(unc%num))

Finally, our required example:

   exp&0.5 (x1 sub x2) add&(exp&2) y1 sub y2
111.803 2.48717

[edit] Java

public class Approx {
private double value;
private double error;
 
public Approx(){this.value = this.error = 0;}
 
public Approx(Approx b){
this.value = b.value;
this.error = b.error;
}
 
public Approx(double value, double error){
this.value = value;
this.error = error;
}
 
public Approx add(Approx b){
value+= b.value;
error = Math.sqrt(error * error + b.error * b.error);
return this;
}
 
public Approx add(double b){
value+= b;
return this;
}
 
public Approx sub(Approx b){
value-= b.value;
error = Math.sqrt(error * error + b.error * b.error);
return this;
}
 
public Approx sub(double b){
value-= b;
return this;
}
 
public Approx mult(Approx b){
double oldVal = value;
value*= b.value;
error = Math.sqrt(value * value * (error*error) / (oldVal*oldVal) +
(b.error*b.error) / (b.value*b.value));
return this;
}
 
public Approx mult(double b){
value*= b;
error = Math.abs(b * error);
return this;
}
 
public Approx div(Approx b){
double oldVal = value;
value/= b.value;
error = Math.sqrt(value * value * (error*error) / (oldVal*oldVal) +
(b.error*b.error) / (b.value*b.value));
return this;
}
 
public Approx div(double b){
value/= b;
error = Math.abs(b * error);
return this;
}
 
public Approx pow(double b){
double oldVal = value;
value = Math.pow(value, b);
error = Math.abs(value * b * (error / oldVal));
return this;
}
 
@Override
public String toString(){return value+"±"+error;}
 
public static void main(String[] args){
Approx x1 = new Approx(100, 1.1);
Approx x2 = new Approx(50, 1.2);
Approx y1 = new Approx(200, 2.2);
Approx y2 = new Approx(100, 2.3);
 
x1.sub(x2).pow(2).add(y1.sub(y2).pow(2)).pow(0.5);
 
System.out.println(x1);
}
}

Output:

111.80339887498948±2.938366893361004

[edit] Mathematica

PlusMinus /: a_ ± σa_ + c_?NumericQ := N[(a + c) ± σa];
PlusMinus /: a_ ± σa_ + b_ ± σb_ := N[(a + b) ± Norm@{σa, σb}];
PlusMinus /: c_?NumericQ (a_ ± σa_) := N[c a ± Abs[c σa]];
PlusMinus /: (a_ ± σa_) (b_ ± σb_) := N[a b ± (a b Norm@{σa/a, σb/b})^2];
PlusMinus /: (a_ ± σa_)^c_?NumericQ := N[a^c ± Abs[a^c σa/a]];
x1 = 100 ± 1.1;
y1 = 50 ± 1.2;
x2 = 200 ± 2.2;
y2 = 100 ± 2.3;
d = Sqrt[(x1 - x2)^2 + (y1 - y2)^2]
Output:
111.803 ± 2.48717

[edit] PARI/GP

This is a work-in-progress.

add(a,b)=if(type(a)==type(b), a+b, if(type(a)=="t_VEC",a+[b,0],b+[a,0]));
sub(a,b)=if(type(a)==type(b), [a[1]-b[1],a[2]+b[2]], if(type(a)=="t_VEC",a-[b,0],[a,0]-b));
mult(a,b)=if(type(a)=="t_VEC", if(type(b)=="t_VEC", [a[1]*b[1], abs(a[1]*b[1])*sqrt(norml2([a[2]/a[1],b[2]/b[1]]))], [b*a[1], abs(b)*a[2]]), [a*b[1], abs(a)*b[2]]);
div(a,b)=if(type(b)!="t_VEC", mult(a,1/b), mult(a,[1/b[1],b[2]/b[1]^2]));
pow(a,b)=[a[1]^b, abs(a[1]^b*b*a[2]/a[1])];
x1=[100,1.1];y1=[50,1.2];x2=[200,2.2];y2=[100,2.3];
pow(add(pow(sub(x1,x2),2),pow(sub(y1,y2),2)),.5)

[edit] Perl

Following code keeps track of covariance between variables. Each variable with error contains its mean value and components of error source from a set of indepentent variables. It's more than what the task requires.

use utf8;
package ErrVar;
use strict;
 
# helper function, apply f to pairs (a, b) from listX and listY
sub zip(&$$) {
my ($f, $x, $y) = @_;
my $l = $#$x;
if ($l < $#$y) { $l = $#$y };
 
my @out;
for (0 .. $l) {
local $a = $x->[$_];
local $b = $y->[$_];
push @out, $f->();
}
\@out
}
 
use overload
'""' => \&_str,
'+' => \&_add,
'-' => \&_sub,
'*' => \&_mul,
'/' => \&_div,
'bool' => \&_bool,
'<=>' => \&_ncmp,
'neg' => \&_neg,
 
'sqrt' => \&_sqrt,
'log' => \&_log,
'exp' => \&_exp,
'**' => \&_pow,
;
 
# make a variable with mean value and a list of coefficient to
# variables providing independent errors
sub make {
my $x = shift;
bless [$x, [@{+shift}]]
}
 
sub _str { sprintf "%g±%.3g", $_[0][0], sigma($_[0]) }
 
# mean value of the var, or just the input if it's not of this class
sub mean {
my $x = shift;
ref($x) && $x->isa(__PACKAGE__) ? $x->[0] : $x
}
 
# return variance index array
sub vlist {
my $x = shift;
ref($x) && $x->isa(__PACKAGE__) ? $x->[1] : [];
}
 
sub variance {
my $x = shift;
return 0 unless ref($x) and $x->isa(__PACKAGE__);
my $s;
$s += $_ * $_ for (@{$x->[1]});
$s
}
 
sub covariance {
my ($x, $y) = @_;
return 0 unless ref($x) && $x->isa(__PACKAGE__);
return 0 unless ref($y) && $y->isa(__PACKAGE__);
 
my $s;
zip { $s += $a * $b } vlist($x), vlist($y);
$s
}
 
sub sigma { sqrt variance(shift) }
 
# to determine if a var is probably zero. we use 1σ here
sub _bool {
my $x = shift;
return abs(mean($x)) > sigma($x);
}
 
sub _ncmp {
my $x = shift() - shift() or return 0;
return mean($x) > 0 ? 1 : -1;
}
 
sub _neg {
my $x = shift;
bless [ -mean($x), [map(-$_, @{vlist($x)}) ] ];
}
 
sub _add {
my ($x, $y) = @_;
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
bless [$x0 + $y0, zip {$a + $b} $xv, $yv];
}
 
sub _sub {
my ($x, $y, $swap) = @_;
if ($swap) { ($x, $y) = ($y, $x) }
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
bless [$x0 - $y0, zip {$a - $b} $xv, $yv];
}
 
sub _mul {
my ($x, $y) = @_;
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
 
$xv = [ map($y0 * $_, @$xv) ];
$yv = [ map($x0 * $_, @$yv) ];
 
bless [$x0 * $y0, zip {$a + $b} $xv, $yv];
}
 
sub _div {
my ($x, $y, $swap) = @_;
if ($swap) { ($x, $y) = ($y, $x) }
 
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
 
$xv = [ map($_/$y0, @$xv) ];
$yv = [ map($x0 * $_/$y0/$y0, @$yv) ];
 
bless [$x0 / $y0, zip {$a + $b} $xv, $yv];
}
 
sub _sqrt {
my $x = shift;
my $x0 = mean($x);
my $xv = vlist($x);
$x0 = sqrt($x0);
$xv = [ map($_ / 2 / $x0, @$xv) ];
bless [$x0, $xv]
}
 
sub _pow {
my ($x, $y, $swap) = @_;
if ($swap) { ($x, $y) = ($y, $x) }
if ($x < 0) {
if (int($y) != $y || ($y & 1)) {
die "Can't take pow of negative number $x";
}
$x = -$x;
}
exp($y * log $x)
}
 
sub _exp {
my $x = shift;
my $x0 = exp(mean($x));
my $xv = vlist($x);
bless [ $x0, [map($x0 * $_, @$xv) ] ]
}
 
sub _log {
my $x = shift;
my $x0 = mean($x);
my $xv = vlist($x);
bless [ log($x0), [ map($_ / $x0, @$xv) ] ]
}
 
"If this package were to be in its own file, you need some truth value to end it like this.";
 
package main;
 
sub e { ErrVar::make @_ };
 
# x1 is of mean value 100, containing error 1.1 from source 1, etc.
# all error sources are independent.
my $x1 = e 100, [1.1, 0, 0, 0 ];
my $x2 = e 200, [0, 2.2, 0, 0 ];
my $y1 = e 50, [0, 0, 1.2, 0 ];
my $y2 = e 100, [0, 0, 0, 2.3];
 
my $z1 = sqrt(($x1 - $x2) ** 2 + ($y1 - $y2) ** 2);
print "distance: $z1\n\n";
 
# this is not for task requirement
my $a = $x1 + $x2;
my $b = $y1 - 2 * $x2;
print "covariance between $a and $b: ", $a->covariance($b), "\n";
output
distance: 111.803±2.49
 
covariance between 300±2.46 and -350±4.56: -9.68

[edit] Perl 6

Translation of: Perl
# cache of independent sources so we can make them all the same length.
# (Because Perl 6 does not yet have a longest-zip metaoperator.)
my @INDEP;
 
class Approx does Numeric {
has Real $.x; # The mean.
has $.c; # The components of error.
 
multi method Str { sprintf "%g±%.3g", $!x, $.σ }
multi method Bool { abs($!x) > $.σ }
 
method variance { [+] @.c X** 2 }
method σ { sqrt self.variance }
}
 
multi approx($x,$c) { Approx.new: :$x, :$c }
multi approx($x) { Approx.new: :$x, :c[0 xx +@INDEP] }
 
# Each ± gets its own source slot.
multi infix:<±>($a, $b) {
.push: 0 for @INDEP; # lengthen older component lists
my $c = [ 0 xx @INDEP, $b ];
@INDEP.push: $c; # add new component list
 
approx $a, $c;
}
 
multi prefix:<->(Approx $a) { approx -$a.x, [$a.c.map: -*] }
 
multi infix:<+>($a, Approx $b) { approx($a) + $b }
multi infix:<+>(Approx $a, $b) { $a + approx($b) }
multi infix:<+>(Approx $a, Approx $b) { approx $a.x + $b.x, [$a.c Z+ $b.c] }
 
multi infix:<->($a, Approx $b) { approx($a) - $b }
multi infix:<->(Approx $a, $b) { $a - approx($b) }
multi infix:<->(Approx $a, Approx $b) { approx $a.x - $b.x, [$a.c Z- $b.c] }
 
multi covariance(Real $a, Real $b) { 0 }
multi covariance(Approx $a, Approx $b) { [+] $a.c Z* $b.c }
 
multi infix:«<=>»(Approx $a, Approx $b) { $a.x <=> $b.x }
multi infix:<cmp>(Approx $a, Approx $b) { $a.x <=> $b.x }
 
multi infix:<*>($a, Approx $b) { approx($a) * $b }
multi infix:<*>(Approx $a, $b) { $a * approx($b) }
multi infix:<*>(Approx $a, Approx $b) {
approx $a.x * $b.x,
[$a.c.map({$b.x * $_}) Z+ $b.c.map({$a.x * $_})];
}
 
multi infix:</>($a, Approx $b) { approx($a) / $b }
multi infix:</>(Approx $a, $b) { $a / approx($b) }
multi infix:</>(Approx $a, Approx $b) {
approx $a.x / $b.x,
[ $a.c.map({ $_ / $b.x }) Z+ $b.c.map({ $a.x * $_ / $b.x / $b.x }) ];
}
 
multi sqrt(Approx $a) {
my $x = sqrt($a.x);
approx $x, [ $a.c.map: { $_ / 2 / $x } ];
}
 
multi infix:<**>(Approx $a, Real $b) { $a ** approx($b) }
multi infix:<**>(Approx $a is copy, Approx $b) {
my $ax = $a.x;
my $bx = $b.x;
my $fbx = floor $b.x;
if $ax < 0 {
if $fbx != $bx or $fbx +& 1 {
die "Can't take power of negative number $ax";
}
$a = -$a;
}
exp($b * log $a);
}
 
multi exp(Approx $a) {
my $x = exp($a.x);
approx $x, [ $a.c.map: { $x * $_ } ];
}
 
multi log(Approx $a) {
my $x0 = $a.x;
approx log($x0), [ $a.c.map: { $_ / $x0 }];
}
 
# Each ± sets up an independent source component.
my $x1 = 100 ± 1.1;
my $x2 = 200 ± 2.2;
my $y1 = 50 ± 1.2;
my $y2 = 100 ± 2.3;
 
# The standard task.
my $z1 = sqrt(($x1 - $x2) ** 2 + ($y1 - $y2) ** 2);
say "distance: $z1\n";
 
# Just showing off.
my $a = $x1 + $x2;
my $b = $y1 - 2 * $x2;
say "covariance between $a and $b: ", covariance($a,$b);
Output:
distance: 111.803±2.49

covariance between 300±2.46 and -350±4.56: -9.68

[edit] PicoLisp

For this task, we overload the built-in arithmetic functions. If the arguments are cons pairs, they are assumed to hold the fixpoint number in the CAR, and the uncertainty's square in the CDR. Otherwise normal numbers are handled as usual.

The overloaded +, -, * and / operators look a bit complicated, because they must handle an arbitrary number of arguments to be compatible with the standard operators.

(scl 12)
(load "@lib/math.l")
 
# Overload arithmetic operators +, -, *, / and **
(redef + @
(let R (next)
(while (args)
(let N (next)
(setq R
(if2 (atom R) (atom N)
(+ R N) # c + c
(cons (+ R (car N)) (cdr N)) # c + a
(cons (+ (car R) N) (cdr R)) # a + c
(cons # a + b
(+ (car R) (car N))
(+ (cdr R) (cdr N)) ) ) ) ) )
R ) )
 
(redef - @
(let R (next)
(ifn (args)
(- R)
(while (args)
(let N (next)
(setq R
(if2 (atom R) (atom N)
(- R N) # c - c
(cons (- R (car N)) (cdr N)) # c - a
(cons (- (car R) N) (cdr R)) # a - c
(cons # a - b
(- (car R) (car N))
(+ (cdr R) (cdr N)) ) ) ) ) )
R ) ) )
 
(redef * @
(let R (next)
(while (args)
(let N (next)
(setq R
(if2 (atom R) (atom N)
(* R N) # c * c
(cons # c * a
(*/ R (car N) 1.0)
(mul2div2 (cdr N) R 1.0) )
(cons # a * c
(*/ (car R) N 1.0)
(mul2div2 (cdr R) N 1.0) )
(uncMul (*/ (car R) (car N) 1.0) R N) ) ) ) ) # a * b
R ) )
 
(redef / @
(let R (next)
(while (args)
(let N (next)
(setq R
(if2 (atom R) (atom N)
(/ R N) # c / c
(cons # c / a
(*/ R 1.0 (car N))
(mul2div2 (cdr N) R 1.0) )
(cons # a / c
(*/ (car R) 1.0 N)
(mul2div2 (cdr R) N 1.0) )
(uncMul (*/ (car R) 1.0 (car N)) R N) ) ) ) ) # a / b
R ) )
 
(redef ** (A C)
(if (atom A)
(** A C)
(let F (pow (car A) C)
(cons F
(mul2div2 (cdr A) (*/ F C (car A)) 1.0) ) ) ) )
 
# Utilities
(de mul2div2 (A B C)
(*/ A B B (* C C)) )
 
(de uncMul (F R N)
(cons F
(mul2div2
(+
(mul2div2 (cdr R) 1.0 (car R))
(mul2div2 (cdr N) 1.0 (car N)) )
F
1.0 ) ) )
 
# I/O conversion
(de unc (N U)
(if U
(cons N (*/ U U 1.0))
(pack
(round (car N) 10)
" ± "
(round (sqrt (cdr N) 1.0) 8) ) ) )

Test:

(de distance (X1 Y1 X2 Y2)
(**
(+ (** (- X1 X2) 2.0) (** (- Y1 Y2) 2.0))
0.5 ) )
 
(prinl "Distance: "
(unc
(distance
(unc 100. 1.1)
(unc 50. 1.2)
(unc 200. 2.2)
(unc 100. 2.3) ) ) )

Output:

Distance: 111.8033988750 ± 2.48716706

[edit] Python

from collections import namedtuple
import math
 
class I(namedtuple('Imprecise', 'value, delta')):
'Imprecise type: I(value=0.0, delta=0.0)'
 
__slots__ = ()
 
def __new__(_cls, value=0.0, delta=0.0):
'Defaults to 0.0 ± delta'
return super().__new__(_cls, float(value), abs(float(delta)))
 
def reciprocal(self):
return I(1. / self.value, self.delta / (self.value**2))
 
def __str__(self):
'Shorter form of Imprecise as string'
return 'I(%g, %g)' % self
 
def __neg__(self):
return I(-self.value, self.delta)
 
def __add__(self, other):
if type(other) == I:
return I( self.value + other.value, (self.delta**2 + other.delta**2)**0.5 )
try:
c = float(other)
except:
return NotImplemented
return I(self.value + c, self.delta)
 
def __sub__(self, other):
return self + (-other)
 
def __radd__(self, other):
return I.__add__(self, other)
 
def __mul__(self, other):
if type(other) == I:
#if id(self) == id(other):
# return self ** 2
a1,b1 = self
a2,b2 = other
f = a1 * a2
return I( f, f * ( (b1 / a1)**2 + (b2 / a2)**2 )**0.5 )
try:
c = float(other)
except:
return NotImplemented
return I(self.value * c, self.delta * c)
 
def __pow__(self, other):
if type(other) == I:
return NotImplemented
try:
c = float(other)
except:
return NotImplemented
f = self.value ** c
return I(f, f * c * (self.delta / self.value))
 
def __rmul__(self, other):
return I.__mul__(self, other)
 
def __truediv__(self, other):
if type(other) == I:
return self.__mul__(other.reciprocal())
try:
c = float(other)
except:
return NotImplemented
return I(self.value / c, self.delta / c)
 
def __rtruediv__(self, other):
return other * self.reciprocal()
 
__div__, __rdiv__ = __truediv__, __rtruediv__
 
Imprecise = I
 
def distance(p1, p2):
x1, y1 = p1
x2, y2 = p2
return ((x1 - x2)**2 + (y1 - y2)**2)**0.5
 
x1 = I(100, 1.1)
x2 = I(200, 2.2)
y1 = I( 50, 1.2)
y2 = I(100, 2.3)
 
p1, p2 = (x1, y1), (x2, y2)
print("Distance between points\n p1: %s\n and p2: %s\n = %r" % (
p1, p2, distance(p1, p2)))
Sample output
Distance between points
  p1: (I(value=100.0, delta=1.1), I(value=50.0, delta=1.2))
  and p2: (I(value=200.0, delta=2.2), I(value=100.0, delta=2.3))
  = I(value=111.80339887498948, delta=2.4871670631463423)

[edit] Racket

Translation of: Mathematica
#lang racket
 
(struct ± (x dx) #:transparent
#:methods gen:custom-write
[(define (write-proc a port mode) (display (±->string a) port))])
 
(define/match (±+ a [b 0])
[((± x dx) (± y dy)) (± (+ x y) (norm dx dy))]
[((± x dx) c) (± (+ x c) dx)]
[(_ (± y dy)) (±+ b a)])
 
(define/match (±* a b)
[((± x dx) (± y dy)) (± (* x y) (* x y (norm (/ dx x) (/ dy y))))]
[((± x dx) c) (± (* x c) (abs (* c dx)))]
[(_ (± y dy)) (±* b a)])
 
(define/match (±- a [b #f])
[(a #f) (±* -1 a)]
[(a b) (±+ a (±* -1 b))])
 
(define/match (±/ a b)
[((± x dx) (± y dy)) (± (/ x y) (/ x y (norm (/ dx x) (/ dy y))))]
[((± _ _) c) (±* a (/ 1 c))])
 
(define/match (±expt a c)
[((± x dx) c) (± (expt x c) (abs (* (expt x c) (/ dx x))))])
 
(define/match (norm a b)
[((± x dx) (± y dy)) (±expt (±+ (±expt a 2) (±expt b 2)) 0.5)]
[(x y) (sqrt (+ (sqr x) (sqr y)))])
 
(define/match (±->string x [places 3])
[((± x dx) p) (string-join (map (λ (s) (real->decimal-string s p))
(list x dx))" ± ")])
 
;; Test
;;
(define x1 (± 100 1.1))
(define y1 (± 50 1.2))
(define x2 (± 200 2.2))
(define y2 (± 100 2.3))
(norm (±- x1 x2) (±- y1 y2))
Output:
111.803 ± 2.487

[edit] Ruby

class NumberWithUncertainty
def initialize(number, error)
@num = number
@err = error.abs
end
attr_reader :num, :err
 
def +(other)
if other.kind_of?(self.class)
self.class.new(num + other.num, Math::hypot(err, other.err))
else
self.class.new(num + other, err)
end
end
 
def -(other)
if other.kind_of?(self.class)
self.class.new(num - other.num, Math::hypot(err, other.err))
else
self.class.new(num - other, err)
end
end
 
def *(other)
if other.kind_of?(self.class)
prod = num * other.num
e = Math::hypot((prod * err / num), (prod * other.err / other.num))
self.class.new(prod, e)
else
self.class.new(num * other, (err * other).abs)
end
end
 
def /(other)
if other.kind_of?(self.class)
quo = num / other.num
e = Math::hypot((quo * err / num), (quo * other.err / other.num))
self.class.new(quo, e)
else
self.class.new(num / other, (err * other).abs)
end
end
 
def **(exponent)
Float(exponent) rescue raise ArgumentError, "not an number: #{exponent}"
prod = num ** exponent
self.class.new(prod, (prod * exponent * err / num).abs)
end
 
def sqrt
self ** 0.5
end
 
def to_s
"#{num} \u00b1 #{err}"
end
end
 
x1 = NumberWithUncertainty.new(100, 1.1)
y1 = NumberWithUncertainty.new( 50, 1.2)
x2 = NumberWithUncertainty.new(200, 2.2)
y2 = NumberWithUncertainty.new(100, 2.3)
 
puts ((x1 - x2) ** 2 + (y1 - y2) ** 2).sqrt

outputs

111.803398874989 ± 2.48716706314634

[edit] Tcl

Works with: Tcl version 8.6

Firstly, some support code for doing RAII-like things, evolved from code in the quaternion solution:

package require Tcl 8.6
oo::class create RAII-support {
constructor {} {
upvar 1 { end } end
lappend end [self]
trace add variable end unset [namespace code {my DieNicely}]
}
destructor {
catch {
upvar 1 { end } end
trace remove variable end unset [namespace code {my DieNicely}]
}
}
method return {{level 1}} {
incr level
upvar 1 { end } end
upvar $level { end } parent
trace remove variable end unset [namespace code {my DieNicely}]
lappend parent [self]
trace add variable parent unset [namespace code {my DieNicely}]
return -level $level [self]
}
# Swallows arguments
method DieNicely args {tailcall my destroy}
}
oo::class create RAII-class {
superclass oo::class
method return args {
[my new {*}$args] return 2
}
method unknown {m args} {
if {[string is double -strict $m]} {
return [tailcall my new $m {*}$args]
}
next $m {*}$args
}
unexport create unknown
self method create args {
set c [next {*}$args]
oo::define $c superclass {*}[info class superclass $c] RAII-support
return $c
}
}
# Makes a convenient scope for limiting RAII lifetimes
proc scope {script} {
foreach v [info global] {
if {[array exists ::$v] || [string match { * } $v]} continue
lappend vars $v
lappend vals [set ::$v]
}
tailcall apply [list $vars [list \
try $script on ok msg {$msg return}
] [uplevel 1 {namespace current}]] {*}$vals
}

The implementation of the number+error class itself:

RAII-class create Err {
variable N E
constructor {number {error 0.0}} {
next
namespace import ::tcl::mathfunc::* ::tcl::mathop::*
variable N $number E [abs $error]
}
method p {} {
return "$N \u00b1 $E"
}
 
method n {} { return $N }
method e {} { return $E }
 
method + e {
if {[info object isa object $e]} {
Err return [+ $N [$e n]] [hypot $E [$e e]]
} else {
Err return [+ $N $e] $E
}
}
method - e {
if {[info object isa object $e]} {
Err return [- $N [$e n]] [hypot $E [$e e]]
} else {
Err return [- $N $e] $E
}
}
method * e {
if {[info object isa object $e]} {
set f [* $n [$E n]]
Err return $f [expr {hypot($E*$f/$N, [$e e]*$f/[$e n])}]
} else {
Err return [* $N $e] [abs [* $E $e]]
}
}
method / e {
if {[info object isa object $e]} {
set f [/ $n [$E n]]
Err return $f [expr {hypot($E*$f/$N, [$e e]*$f/[$e n])}]
} else {
Err return [/ $N $e] [abs [/ $E $e]]
}
}
method ** c {
set f [** $N $c]
Err return $f [abs [* $f $c [/ $E $N]]]
}
 
export + - * / **
}

Demonstrating:

set x1 [Err 100 1.1]
set x2 [Err 200 2.2]
set y1 [Err 50 1.2]
set y2 [Err 100 2.3]
# Evaluate in a local context to clean up intermediate objects
set d [scope {
[[[$x1 - $x2] ** 2] + [[$y1 - $y2] ** 2]] ** 0.5
}]
puts "d = [$d p]"

Output:

d = 111.80339887498948 ± 2.4871670631463423
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox