{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--loading configuration for package \"FourTiTwo\" from file /home/hegland/.Macaulay2/init-FourTiTwo.m2\n",
"--loading configuration for package \"Topcom\" from file /home/hegland/.Macaulay2/init-Topcom.m2\n"
]
}
],
"source": [
"restart"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rounding\n",
"\n",
"One of the challenges of rational arithmetic is that the numerators and denominators can grow strongly during computations. One way to control this growth uses approximate computations and a rounding function. The rounding function is a real function and its range is the set $\\mathbb{F}_t$ of floating point numbers, a subset of the ring of dyadic fractions (aka dyadic rationals).\n",
"\n",
"Dyadic fractions are the extension $\\mathbb{Z}[1/2]$ of the ring of integers by the fraction $\\frac{1}{2}$. Any dyadic fraction can be (non-uniquely) represented by a mantissa $m\\in\\mathbb{Z}$ and an exponent $e\\in\\mathbb{N}$ as \n",
" $$x = m\\, 2^{-e}.$$\n",
"The subset of dyadic fractions with $e=0$ is the set of integers. One thus has\n",
" $$\\mathbb{Z} \\subset \\mathbb{Z}[1/2] \\subset \\mathbb{Q}.$$\n",
"Like the rationals the dyadic fractions are dense in the set of real numbers. While they are not a field, they require around half the storage required for a rational number as the exponent typically is much smaller and thus has much less decimal digits than the denominator of a fraction.\n",
"\n",
"The set of floating point numbers $\\mathbb{F}_t$ is obtained by limiting the range of the mantissa:\n",
" $$\\mathbb{F}_t = \\{0\\}\\cup\\{m\\, 2^e \\mid 2^{t-1} \\leq \\lvert m \\rvert < 2^t,\\; e,m\\in\\mathbb{Z}\\}.$$\n",
" This representation of the floating point numbers is called \"normalised\". It is unique.\n",
"The set of floating point numbers is not a ring and has one accumulation point 0. The set of floating point numbers contains a subset of the set of integers (not all of them) and is a subset of the set of dyadic fractions.\n",
"\n",
"* The set $\\mathbb{F}_t$ contains all integers between $-2^t$ and $2^t$.\n",
"\n",
" Proof: Choosing $e=0$ on sees that it contains the integers between $2^{t-1}$ and $2^t$ and $e=-1...-(t-1)$ revealls that all the positive integers less than $2^{t-1}$ are also contained. Finally the set is symmetric and contains 0.\n",
" \n",
"* Let $\\mathbb{Z}_{2^t+1}=\\{-2^{t-1},\\ldots, 2^{t-1}\\}$. Then\n",
"\n",
" $$\\mathbb{F}_t = \\{0\\}\\cup\\bigcup_{e\\in\\mathbb{Z}} 2^e\\mathbb{Z}_{2^{t+1}+1}.$$\n",
" This leads to representations of the floating point numbers which are not unique.\n",
" \n",
"* Consequently $2\\mathbb{F}_t = \\mathbb{F}_t$.\n",
"\n",
"The rounding function $\\rho_t : \\mathbb{R} \\rightarrow \\mathbb{F}_t$ is based on a rounding function $\\rho_0 : \\mathbb{R} \\rightarrow \\mathbb{Z}$ which is provided by Macaulay. This function returns the closest integer to any given real argument $x$. In particular, this function satisfies $\\rho_0(x)=0$ if $x\\in[0,1/2)$ and $\\rho_0(x)=1$ if $x\\in[1/2,1)$.\n",
"\n",
"Now let $x=y\\, 2^e$ for positive $x\\in\\mathbb{R}$ be such that\n",
"$y\\in[2^{t-1},2^t)$. (Like the normalised floating point representation, this representation of a positive real number is unique.) Then we define $\\rho_t(x)$ to be\n",
" $$\\rho_t(x) = \\rho_0(y)\\, 2^e.$$\n",
"This representation is not normalised if $\\rho_0(y)=2^t$. Finally we set $\\rho_t(0)=0$ and $\\rho_t(-x) = -\\rho_t(x)$.\n",
"\n",
"Then $\\rho_t$ has the following properties:\n",
"\n",
"* $\\rho_t(2^s x) = 2^s \\rho_t(x)$ for any $s\\in\\mathbb{Z}$.\n",
"\n",
" For example one has $\\rho_3(7.1) = 7$ as $7.1\\in [4,8)$ thus $7.1=7.1*2^0$. Then $3.55= 7.1*2^{-1}$ and thus $\\rho_3(7.1/2)=7/2$. The rounding function $\\rho_0$ does not satisfy this property in general. One has $\\rho_0(7.1) = 7 \\neq 2*\\rho_0(3.55) = 8$.\n",
" \n",
" The property follows directly from the uniqueness of the representation of real positive numbers as $x=z\\,2^e$ with $z\\in[2^{t-1},2^t)$.\n",
" \n",
"* For all $x\\neq 0$ one has\n",
" $$\\frac{\\lvert \\rho_t(x)-x\\rvert}{\\lvert x \\rvert} \\leq 2^{-t}.$$\n",
" \n",
" Proof: As in the definition let $x=y 2^e$ where $2^{t-1}\\leq \\lvert y \\rvert < 2^t$.\n",
" Substituting the definition of $\\rho_t(x)$ one then has \n",
" $$\\frac{\\lvert \\rho_t(x)-x\\rvert}{\\lvert x \\rvert} =\n",
" \\frac{\\lvert \\rho_0(y)-y\\rvert}{\\lvert y \\rvert}$$\n",
" As the rounding error of $\\rho_0$ satisfies\n",
" $$|\\rho_0(y) - y| \\leq {1 \\over 2}$$\n",
" and as $\\lvert y \\rvert \\geq 2^{t-1}$ one gets the claimed\n",
" inequality. QED\n",
" \n",
" Thus one has a bound for the relative error (the \"number of digits\") in contrast to $\\rho_0$ which has a bounded absolute\n",
" error of 1/2.\n",
" \n",
"We can now rewrite the error bound as\n",
" $$\\rho_t(x) - x = \\delta x$$\n",
"or\n",
" $$\\rho_t(x) = (1+\\delta)x$$\n",
"for some real number $\\delta$ with $\\lvert\\delta\\rvert \\leq 2^{-t}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### An algebraic model for uncertainty\n",
"\n",
"* the following code is for illustrative purposes and might\n",
" contain (coding) errors"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"o2 = Rd\n",
"\n",
"o2 : PolynomialRing\n"
]
}
],
"source": [
"-- choose the ring for your computations\n",
"R = QQ -- rational numbers\n",
"Rd = R[delta_1..delta_4] -- delta_i parameters with uncertain value"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"o6 = Ru\n",
"\n",
"o6 : QuotientRing\n"
]
}
],
"source": [
"-- the ring for computations with uncertainty\n",
"-- here we consider affine polynomials of degree one\n",
"nu = numgens Rd\n",
"Lu = flatten(for i from 0 to nu-1 list \n",
" (for j from i to nu-1 list Rd_i*Rd_j))\n",
"Iu = ideal(Lu)\n",
"Ru = Rd/Iu"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" 3 10 4\n",
"o7 = -delta + delta + --delta + -delta + 10\n",
" 2 1 2 7 3 9 4\n",
"\n",
"o7 : Rd\n"
]
}
],
"source": [
"-- a random uncertain expression\n",
"x = random(R) + sum(for i from 0 to nu-1 list Rd_i*random(R))"
]
},
{
"cell_type": "code",
"execution_count": 140,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"expected value = 10\n",
"\n",
"variance factor = 5.48835\n",
"\n",
"max error factor = 1.5\n"
]
}
],
"source": [
"-- the expected value\n",
"mu = coefficient(1_Rd,x)\n",
"err = x - mu\n",
"<< \"expected value = \" << promote(mu,RR) << endl;\n",
"\n",
"-- the variance (uncertain variables are N(0,eps^2))\n",
"var = sum(for i from 0 to nu-1 list (coefficient(Rd_i,x))^2)\n",
"<< \"variance factor = \" << promote(var,RR) << endl;\n",
"\n",
"-- the maximal error\n",
"--errmax = max(for i from 0 to nu-1 list abs(C_0_i))\n",
"errmax = max(for i from 0 to nu-1 list abs(coefficient(Rd_i,x)))\n",
"<< \"max error factor = \" << promote(errmax,RR) << endl;"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"o122 = -1.85037170770859e-17\n",
"\n",
"o122 : RR (of precision 53)\n"
]
}
],
"source": [
"-- the rounding function which creates uncertainty\n",
"\n",
"rho = (x,t,DF) -> (\n",
" if x == 0 then return 0_DF\n",
" else if x > 0 then (\n",
" m = x; f=1_DF; \n",
" while m < 2^(t-1) do (m=2*m; f=h*f);\n",
" while m >= 2^t do (m=m/2; f = 2*f);\n",
" return round(m)*f)\n",
" else return -rho(-x,t,h))\n",
"\n",
"FR(rho(1/3,53,DF))-1/3 -- same as for floating point RR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### lifting functions defined on $\\mathbb{D}$ to $\\mathbb{F}$\n",
"\n",
"* simple formula for unary functions\n",
" $$f_\\mathbb{F}(x) = \\rho(E(f_\\mathbb{D}(x))$$\n",
"* this is an approximation\n",
"* the same for binary functions\n",
" $$f_\\mathbb{F}(x,y) = \\rho(E(f_\\mathbb{D}(x,y)))$$\n",
"* here $E$ is the embedding of $\\mathbb{D}$ into $\\mathbb{Q}$,\n",
" i.e., the function FQ from above\n",
"* one could also implement the rounding function on $\\mathbb{D}$\n",
" instead of $\\mathbb{Q}$"
]
},
{
"cell_type": "code",
"execution_count": 131,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"error in addition of 1/3+1/4 = .03125\n",
"\n",
"sum of errors in approx of 1/3 and 1/4 =.0104167\n"
]
}
],
"source": [
"t = 4\n",
"UF = f -> (x -> rho(FQ(f(x)),t,DF)) -- unary functions f\n",
"BF = f -> ((x,y) -> rho(FQ(f(x,y)),t,DF))-- binary functions f\n",
"\n",
"x = rho(1/3, t, DF) -- rho incurs an error\n",
"y = rho(1/4, t, DF) -- rho incurse an error\n",
"errplus = FR((BF(plus))(x,y) - (x+y)) -- BF(plus) incurs an error\n",
"<< \"error in addition of 1/3+1/4 = \"<< errplus << endl;\n",
"<< \"sum of errors in approx of 1/3 and 1/4 =\"<< FR(x+y)-(1/3+1/4) <