restart
-- choose the ring for your computations
R = QQ -- rational numbers
Rd = R[delta_1..delta_4] -- delta_i parameters with uncertain value
-- the ring for computations with uncertainty
-- here we consider affine polynomials of degree one
nu = numgens Rd
Lu = flatten(for i from 0 to nu-1 list
(for j from i to nu-1 list Rd_i*Rd_j))
Iu = ideal(Lu)
Ru = Rd/Iu
-- a random uncertain expression
x = random(R) + sum(for i from 0 to nu-1 list Rd_i*random(R))
-- the expected value
mu = coefficient(1_Rd,x)
err = x - mu
<< "expected value = " << promote(mu,RR) << endl;
-- the variance (uncertain variables are N(0,eps^2))
var = sum(for i from 0 to nu-1 list (coefficient(Rd_i,x))^2)
<< "variance factor = " << promote(var,RR) << endl;
-- the maximal error
--errmax = max(for i from 0 to nu-1 list abs(C_0_i))
errmax = max(for i from 0 to nu-1 list abs(coefficient(Rd_i,x)))
<< "max error factor = " << promote(errmax,RR) << endl;
-- the rounding function which creates uncertainty
rho = (x,t,DF) -> (
if x == 0 then return 0_DF
else if x > 0 then (
m = x; f=1_DF;
while m < 2^(t-1) do (m=2*m; f=h*f);
while m >= 2^t do (m=m/2; f = 2*f);
return round(m)*f)
else return -rho(-x,t,h))
FR(rho(1/3,53,DF))-1/3 -- same as for floating point RR
t = 4
UF = f -> (x -> rho(FQ(f(x)),t,DF)) -- unary functions f
BF = f -> ((x,y) -> rho(FQ(f(x,y)),t,DF))-- binary functions f
x = rho(1/3, t, DF) -- rho incurs an error
y = rho(1/4, t, DF) -- rho incurse an error
errplus = FR((BF(plus))(x,y) - (x+y)) -- BF(plus) incurs an error
<< "error in addition of 1/3+1/4 = "<< errplus << endl;
<< "sum of errors in approx of 1/3 and 1/4 ="<< FR(x+y)-(1/3+1/4) <