Recall the grade-school approximation to the number π.

@show p = 22/7;

@show float(π);

acc = abs(p-π)
println("absolute accuracy = $acc")
println("relative accuracy = $(acc/π)")

println("Number of accurate digits = $(-log10(acc/π))")

This last value could be rounded down by using floor.

In Julia, 1 and 1.0 are different values, because they have different types:

@show typeof(1);
@show typeof(1.0);

The standard choice for floating-point values is Float64, which is double precision using 64 binary bits. We can see all the bits by using bitstring.

bitstring(1.0)

:::{card}

Square brackets concatenate the contained values into vectors.

[bitstring(1.0), bitstring(-1.0)]

The next 11 bits determine the exponent (scaling) of the number, and so on.

[bitstring(1.0), bitstring(2.0)]

The sign bit, exponent, and significand in (1.1.1) are all directly accessible.

x = 3.14
@show sign(x), exponent(x), significand(x);

x = x / 8
@show sign(x), exponent(x), significand(x);

eps()

Because double precision allocates 52 bits to the significand, the default value of machine epsilon is 2^-52.

log2(eps())

The spacing between adjacent floating-point values is proportional to the magnitude of the value itself. This is how relative precision is kept roughly constant throughout the range of values. You can get the adjusted spacing by calling eps with a value.

eps(1.618)

eps(161.8)

nextfloat(161.8)

ans - 161.8

A common mistake is to think that $\epsilon_\text{mach}$ is the smallest floating-point number. It’s only the smallest relative to 1. The correct perspective is that the scaling of values is limited by the exponent, not the mantissa. The actual range of positive values in double precision is

@show floatmin(),floatmax();

For the most part you can mix integers and floating-point values and get what you expect.

1/7

37.3 + 1

2^(-4)

There are some exceptions. A floating-point value can’t be used as an index into an array, for example, even if it is numerically equal to an integer. In such cases you use Int to convert it.

@show 5.0, Int(5.0);

If you try to convert a noninteger floating-point value into an integer you get an InexactValue error. This occurs whenever you try to force a type conversion that doesn’t make clear sense.

There is no double-precision number between 1 and $1+\epsilon_\text{mach}$. Thus the following difference is zero despite its appearance.

e = eps()/2
(1.0 + e) - 1.0

However, the spacing between floats in $[1/2,1)$ is $\macheps/2$, so both $1-\macheps/2$ and its negative are represented exactly:

1.0 + (e - 1.0)

This is now the expected result. But we have found a rather shocking breakdown of the associative law of addition!

ϵ = 1e-6   # type \epsilon and then press Tab
a,b,c = 1/3,(-2-ϵ)/3,(1+ϵ)/3   # coefficients of p

Here are the roots as computed by the quadratic formula.

d = sqrt(b^2 - 4a*c)
r₁ = (-b - d) / (2a)   # type r\_1 and then press Tab
r₂ = (-b + d) / (2a)
(r₁, r₂)

The relative errors in these values are

@show abs(r₁ - 1) / abs(1);
@show abs(r₂ - (1+ϵ)) / abs(1+ϵ);

The condition number of each root is

Thus, relative error in the data at the level of roundoff can grow in the result to be roughly

eps() / ϵ

This matches the observation pretty well.

Here we show how to use Function 1.3.1 to evaluate a polynomial. It’s not a part of core Julia, so you need to download and install this text’s package once, and load it for each new Julia session. The download is done by the following lines.

import Pkg
Pkg.add("FundamentalsNumericalComputation");

using FundamentalsNumericalComputation

:::{card}

If you are not sure where a particular function is defined, you can run `methods` on the function name to find all its definitions.

Returning to horner, let us define a vector of the coefficients of $p(x)=(x-1)^3=x^3-3x^2+3x-1$, in ascending degree order.

c = [-1, 3, -3, 1]

:::{card}

You must use the module name when a package is loaded by `import`, but when loaded via `using`, some functions may be available with no prefix.

FNC.horner(c, 1.6)

The above is the value of $p(1.6)$.

While the namespace does lead to a little extra typing, a nice side effect of using this paradigm is that if you type FNC. (including the period) and hit the Tab key, you will see a list of all the functions known in that namespace.

The multi-line string at the start of Function 1.3.1 is documentation, which we can access using ?FNC.horner.

a = 1;  b = -(1e6+1e-6);  c = 1;
@show x₁ = (-b + sqrt(b^2 - 4a*c)) / 2a;
@show x₂ = (-b - sqrt(b^2 - 4a*c)) / 2a;

The first value is correct to all stored digits, but the second has fewer than six accurate digits:

error = abs(1e-6-x₂)/1e-6 
@show accurate_digits = -log10(error);

The instability is easily explained. Since $a=c=1$, we treat them as exact numbers. First, we compute the condition numbers with respect to $b$ for each elementary step in finding the “good” root:

Using (1.2.9), the chain rule for condition numbers, the conditioning of the entire chain is the product of the individual steps, so there is essentially no growth of relative error here. However, if we use the quadratic formula for the “bad” root, the next-to-last step becomes $u_4=(-u_3) - b$, and now $\kappa=|u_3|/|u_4|\approx 5\times 10^{11}$. So we can expect to lose 11 digits of accuracy, which is what we observed. The key issue is the subtractive cancellation in this one step.

We repeat the rootfinding experiment of Demo 1.4.1 with an alternative algorithm.

a = 1;  b = -(1e6+1e-6);  c = 1;

First, we find the “good” root using the quadratic formula.

@show x₁ = (-b + sqrt(b^2-4a*c)) / 2a;

Then we use the identity $x_1x_2=\frac{c}{a}$ to compute the smaller root:

@show x₂ = c / (a*x₁);

As you see in this output, Julia often suppresses trailing zeros in a decimal expansion. To be sure we have an accurate result, we compute its relative error.

abs(x₂-1e-6) / 1e-6

using FundamentalsNumericalComputation

Our first step is to construct a polynomial with six known roots.

r = [-2.0,-1,1,1,3,6]
p = fromroots(r)

Now we use a standard numerical method for finding those roots, pretending that we don’t know them already. This corresponds to $\tilde{y}$ in Definition 1.4.1.

r̃ = sort(roots(p))   # type r\tilde and then press Tab

println("Root errors:") 
@. abs(r-r̃) / r

It seems that the forward error is acceptably close to machine epsilon for double precision in all cases except the double root at $x=1$. This is not a surprise, though, given the poor conditioning at such roots.

Let’s consider the backward error. The data in the rootfinding problem is the polynomial coefficients. We can apply fromroots to find the coefficients of the polynomial (that is, the data) whose roots were actually computed by the numerical algorithm. This corresponds to $\tilde{x}$ in Definition 1.4.1.

p̃ = fromroots(r̃)

We find that in a relative sense, these coefficients are very close to those of the original, exact polynomial:

c,c̃ = coeffs(p),coeffs(p̃)
println("Coefficient errors:") 
@. abs(c-c̃) / c

In summary, even though there are some computed roots relatively far from their correct values, they are nevertheless the roots of a polynomial that is very close to the original.