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

 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.
 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)²   +   (y1 - y2)²    
 Print and display both    d    and its error.


References

 A Guide to Error Propagation B. Keeney, 2005.
 Propagation of uncertainty Wikipedia.


Related task

   Quaternion type


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";