Give an exact formula for the smallest positive intiger n that does not belong to $\mathbb{F}$
$\mathbb{F}$ is made up of numbers in the form
$$ x = \pm (m/ \beta^t) \beta^e$$We know the value $\beta^t$ must be in $\mathbb{F}$ because that value can be represented with the significand alone. Assuming $\beta \ge 2$, smallest value not in $\mathbb{F}$ is $\beta^t + 1$, because all values between $\beta^t$ and $\beta^{t + 1}$ are multiplied by two, so odd values are not represented.
What are those values of n for IEEE single and double precision arithmetic?
-Double-precision binary floating-point form has 53 digits of significand stored so $n = 2^{53} +1$ is not stored exactly.
-Double-precision binary floating-point form has 24 digits of significand stored so $n = 2^{24} +1$ is not stored exactly.
In [32]:
ex = 53
m64 = float64((2^ex))
#nbig = big(2^e)
mbig = BigInt((2^ex))
### verify that n = 2^{54} +1 is not represented by float64
(m64 +1) == (mbig+1)
Out[32]:
In [33]:
### verify that n -1 = 2^{54}, n -2 = 2^{54} -1, and n -2 = 2^{54} -2 are represented by float64
(m64+0) == (mbig+0),(m64 -1) == (mbig-1),(m64 -2) == (mbig-2)
Out[33]:
In [34]:
(m64+2) == (mbig+2),(m64 + 4) == (mbig + 4),(m64 +6) == (mbig + 6)
Out[34]:
In [35]:
(m64+1) == (mbig+1),(m64 + 3) == (mbig + 3),(m64 +5) == (mbig + 5)
Out[35]:
In [4]:
b =1
c = logspace(-20,-1,50)
length(c)
x = b - sqrt(b^2 - c)
# multiply by b^2
set_bigfloat_precision(256)
xbig = b - sqrt(b^2 - big(c))
err = float64(abs(x-xbig) ./ abs(xbig))
using PyPlot
loglog(c, err, "bo-")
xlabel(L"c")
ylabel(L"relative error in $x$")
title(L"standard quadratic formula $x = b - \sqrt{b^2 - c}$")
Out[4]:
In [158]:
b =1
c = logspace(-20,-1,50)
length(c)
x = c ./ (b + sqrt(b^2 - c))
set_bigfloat_precision(256)
xbig = b - sqrt(b^2 - big(c))
err = float64(abs(x-xbig) ./ abs(xbig))
loglog(c, err, "bo-")
xlabel(L"c")
ylabel(L"relative error in $x$")
title(L"Alternative calculation of root: $\frac{c}{b + \sqrt{b^2 - c}}$")
Out[158]:
*The last equation represents an iteration for a Newton-like root estimation
Analyze its asymptotic convergence rate: if x is an exact root, write $x_n = x(1 + δn)$ as in class, and solve for δn+1 in terms of δn assuming you are close to the root (δn << 1).
If a single iteration of our Newtonish method takes the form
$$ x_{n+1} = x_n - \frac{2f(x_n) f'(x_n)}{2f'(x_n)^2 - f''(x_n)(f(x_n))} \\ $$Then estimate at a given Newtonish step will have an error associated with it $x_n = x(1 + δ_n)$ and so will the estimate on the next iteration $x_{n+1} = x(1 + δ_{n+1})$. We can plug these in:
$$ x(1 + δ_{n+1}) = x(1 + δ_n) - \frac{2f(x(1 + δ_n)) f'(x(1 + δ_n))}{2f'(x(1 + δ_n))^2 - f''(x(1 + δ_n))(f(x(1 + δ_n)))} \\ $$Then solve for $δ_{n+1}$
$$ δ_{n+1} = δ_n - \frac{2f(x(1 + δ_n)) f'(x(1 + δ_n))}{2xf'(x(1 + δ_n))^2 - xf''(x(1 + δ_n))(f(x(1 + δ_n)))} \\ $$Modify the Julia Newton’s-method notebook from class to implement your method to compute a root of $f(x) = x^3 − 1$. In particular start with x1 = 2, so that your scheme should(!) converge to x = 1, and look at the error $x_n − 1$. Demonstrate that it agrees with your predicted convergence rate from the previous part. [You should use arbitrary precision as in the notebook from class, so that you can watch the convergence rate for many digits. An approximate number of accurate digits is given by − log10(xn − 1).]
In [48]:
function iterate_newton_basic(n,x)
# find root for x^3−1
# n is number of iterations
# x is initial guess
for i = 1:n
x = x - (x^3 - 1)/(3x^2)
end
return x
end
function iterate_newton_like(n,x)
# find root for x^3−1
# n is number of iterations
# x is initial guess
for i = 1:n
x = x - ( (6*x^5) - (6*x^2) )/( (12*x^4) + 6*x)
end
return x
end
Out[48]:
In [51]:
x1 = 2.0 #initial guess
iterations = 15
xn_basic = zeros(iterations)
xn_like = zeros(iterations)
err_basic = zeros(iterations)
err_like = zeros(iterations)
for n = 1:iterations
xn_basic[n] = iterate_newton_basic(n,x1)
xn_like[n] = iterate_newton_like(n,x1)
err_basic[n] = iterate_newton_basic(n,x1) -1
err_like[n] = iterate_newton_like(n,x1) -1
end
plot(1:iterations, err_basic, "bo-",label="basic Newton")
plot(1:iterations, err_like, "ro-",label="Newtonish")
#loglog(1:iterations, xn_like, "ro-")
xlabel("iteration number")
ylabel(L"error in $x$")
title(L"convergence rate for $x^3 -1$")
legend(loc="upper right",fancybox="true")
Out[51]:
Show that algorithm $\widetilde{f}(x)$ for the function $ f(x) = \sum_{i=1}^{n} x_i $: $$ \widetilde{f}(x) = \sum_{i=1}^{n} x_i \prod_{k=i}^{n} (1+\varepsilon_k ) \\ \text{where } \varepsilon_0 = 0 \text{ and } |\varepsilon_k| \le \varepsilon_{machine} $$
Start with an inductive definition of the sum: $s_0 = 0 $ and $s_k = s_{k-1} + x_k $ for $ 0 \le k \le n$
If $ f(x) = s_n $ then $ \tilde{f}(x) = \tilde{s}_n $ where $$ \tilde{s}_k = \tilde{s}_{k-1} \oplus x_k $$ $$ = (\tilde{s}_{k-1} + x_k) (1 + \epsilon_k)$$
Expanding this form:
$ \tilde{f}(x) = \tilde{s}_n = (\tilde{s}_{n-1} + x_n) (1 + \epsilon_n)$
$= (( \tilde{s}_{n-2} \oplus x_{n-1}) + x_n)(1+\epsilon_n)$
$ = ((( \tilde{s}_{n-3} \oplus x_{n-2}) + x_{n-1})(1+ \epsilon_{n-1}) + x_n)(1+\epsilon_n) $
$ = (((( \tilde{s}_{n-4} \oplus x_{n-3}) + x_{n-2})(1 +\epsilon_{n-2}) + x_{n-1})(1+ \epsilon_{n-1}) + x_n)(1+\epsilon_n) $
Continuing in this form, breaking each $\tilde{s}_k$ down to the base case where we can solve for $\tilde{s}_2$:
$ \tilde{s}_2 = \tilde{s}_{1} \oplus x_2 = (x_1(1 + \epsilon_1) + x_2)(1 + \epsilon_2) $
The general form can be rearranged
$ \tilde{f}(x) = ( \tilde{s}_{n-4} \oplus x_{n-3})(1+ \epsilon_{n-2})(1+ \epsilon_{n-1})(1+ \epsilon_{n}) + x_{n-2}(1 + \epsilon_{n-2})(1+ \epsilon_{n-1})(1+ \epsilon_{n}) + x_{n-1}(1+ \epsilon_{n-1})(1 + \epsilon_{n}) + x_n(1+\epsilon_n)$
This is equivent to:
$ = x_1(1+ \epsilon_{1})(1+ \epsilon_{2})\dots (1+ \epsilon_{n}) + x_2(1+ \epsilon_{2})\dots (1+ \epsilon_{n}) + \dots + x_{n-1}(1+ \epsilon_{n-1})(1 + \epsilon_{n}) + x_n(1+\epsilon_n)$
From here we obtain the final form: $$ \widetilde{f}(x) = \sum_{i=1}^{n} x_i \prod_{k=i}^{n} (1+\varepsilon_k ) \\ \text{where } \varepsilon_0 = 0 \text{ and } |\varepsilon_k| \le \varepsilon_{machine} $$
Show that the error can be bounded by $$ | \widetilde{f}(x) - f(x) | \le n \varepsilon_{machine} \sum_{i=1}^{n} |x_i| + O(\varepsilon_{mach}^2) $$
$$ \text{where } \varepsilon_0 = 0 \text{ , } |\varepsilon_k| \le \varepsilon_{machine} \text{ and }\\ \widetilde{f}(x) = \sum_{i=1}^{n} x_i \prod_{k=i}^{n} (1+\varepsilon_k ) \\ $$The term $\prod_{k=i}^{n} (1+\varepsilon_k )$ can be expanded as follows:
$$\prod_{k=i}^{n} (1+\varepsilon_k ) = (1+ \epsilon_i)(1+\epsilon_{i+1})...(1+\epsilon_n) $$Distributing and accounting that $ |\varepsilon_k| \le \varepsilon_{machine} $ $$\prod_{k=i}^{n} (1+\varepsilon_k ) \le 1 + \epsilon_i + \epsilon_{i+1} + ... \epsilon_n + O(\epsilon_{mach}^2) \\ = 1 + \sum_{i}^n \epsilon_i + O(\epsilon_{mach}^2) $$
Letting $\alpha_i = \sum_{i}^n \epsilon_i + O(\epsilon_{mach}^2)$
$$ \widetilde{f}(x) = \sum_{i}^{n} x_i (1 + \alpha_i) \text { and} \\ | \widetilde{f}(x) - f(x) |= |\sum_{i}^{n} x_i (1 + \alpha_i) - \sum_{i}^{n} x_i| \\ = |\sum_{i}^{n} x_i \alpha_i| \\ = | \sum_{i}^{n} x_i (\sum_{i}^n \epsilon_i + O(\epsilon_{mach}^2)) | \\ \le \sum_{i}^{n} |x_i| n \epsilon_{mach} + O(\epsilon_{mach}^2) $$Rearranging yields
$$ | \widetilde{f}(x) - f(x) | \le n \varepsilon_{machine} \sum_{i=1}^{n} |x_i| + O(\varepsilon_{mach}^2) $$The main result for the random walk is that the expected number of equalizations $g_{2m}$ is proportional to the length of the walk $2m$
$$ g_{2m} \sim \sqrt{\frac{2}{\pi}} \sqrt{2m} $$In errors accumulating from sequential summation, the number of additions n is analogous to the length traveled $2m$ during a random walk. This result explains the mean error bound for randomly distributed values of $\varepsilon_{k}$ because $\sqrt{n}$ can be used to represent the average error away from the 0 that accumulates during n additions. This yields the result
$$ | \tilde{f}(x) - f(x) | = O( \sqrt{n} \varepsilon_{machine} \sum_{i=1}^{n} x_i)) $$In randomly distributed values of $\varepsilon_{k}$ positive and negative can cancel out in some cases, yielding an overall lower bound.
For simplicity, assume n is a power of 2 (so that the set of numbers to add divides evenly in two at each stage of the recursion).Show that the worst-case error bound grows logarithmically rather than linearly with n!
$$ \tilde{f}(x) = \left\{\begin{matrix} 0 & \text{ if } n=0 \\ x_1 & \text{ if } n = 1\\ \tilde{f}(x_{1: \lfloor n/2 \rfloor }) \oplus \tilde{f}(x_{ \lfloor n/2 \rfloor +1 :n }) + & \text{ if } n > 1 \end{matrix}\right. $$Assume n is in the form $n=2^m$ where $m = \log_{2}{n}$. In this case $\lfloor n/2 \rfloor = 2^{m-1}$. Therefore for n>1:
$$ \tilde{f}(x) = \tilde{f}(x_{1:2^{m-1}}) \oplus \tilde{f}(x_{2^{m-1}+1 : 2^m}) $$In sequence summation, the $x_1$ entry takes part in n-1 number of additions, the $x_2$ entry takes part in n-2 number of additions, and so on, while the $x_n$ entry is only added once. This is not the case for recursive summation where each entry takes part in the exactly $\log_{2}{n}$ additions. This leads to a similar form as Problem 4:
$$ \tilde{f}(x) \le \sum_{i=1}^{n} x_i \prod_{k=1}^{\log_{2}{n}} (1+\varepsilon_{mach}) $$Similar to part 4 the term $\prod_{k=i}^{\log_{2}{n}} (1+\varepsilon_k )$ can be expanded by distributing and accounting that $ |\varepsilon_k| \le \varepsilon_{machine} $ $$ \prod_{k=i}^{\log_{2}{n}} (1+\varepsilon_k ) \le 1 + \log_{2}{n} \epsilon_{mach} + O(\epsilon_{mach}^2) $$
Letting $\alpha_i = \log_{2}{n} \epsilon_{mach} + O(\epsilon_{mach}^2)$
$$ \widetilde{f}(x) \le \sum_{i}^{n} x_i (1 + \alpha_i) \text { and} \\ | \widetilde{f}(x) - f(x) | \le |\sum_{i}^{n} x_i (1 + \alpha_i) - \sum_{i}^{n} x_i| \\ = |\sum_{i}^{n} x_i \alpha_i| \\ \le \sum_{i}^{n} |x_i| \log_{2}{n} \epsilon_{mach} + O(\epsilon_{mach}^2) $$Rearranging yields
$$ | \widetilde{f}(x) - f(x) | \le \log_{2}{n} * \varepsilon_{machine} \sum_{i=1}^{n} |x_i| + O(\varepsilon_{mach}^2) $$Similar to Part 4c, total accumulated error after the addition should be proportional to the square root of the number of additions performed. For the sequential addition scheme, there were $n$ additions performed so error was proportional to $\sqrt{n}$. In the recursive scheme, each entry takes part in $\log_{2}{n}$ additions, so the overall error is proportional to $\sqrt{\log_{2}{n}}$
$$ | \tilde{f}(x) - f(x) | = O( \sqrt{\log_{2}{n}} * \varepsilon_{machine} \sum_{i=1}^{n} x_i)) $$In the pset 1 Julia notebook, there is a function “div2sum” that computes ˜f(x) =div2sum(x) in single precision by the above algorithm. Modify it to not be horrendously slow via your suggestion in (c), and then plot its errors for random inputs as a function of n with the help of the example code in the Julia notebook (but with a larger range of lengths n). Are your results consistent with your error estimates above?
In [53]:
# Sum x[first:last]. This function works, but is a little slower than we would like.
function div2sum(x, first=1, last=length(x))
n = last - first + 1;
if n < 2
s = zero(eltype(x))
for i = first:last
s += x[i]
end
return s
else
mid = div(first + last, 2) # find middle as (first+last)/2, rounding down
return div2sum(x, first, mid) + div2sum(x, mid+1, last)
end
end
# make function faster
function div2sum_performance(x, first=1, last=length(x))
n = last - first + 1;
while n > 1
xsum = zeros(eltype(x),int(n/2))
for i = 1:int(n/2)
odd_val = 2i-1
xsum[i] = x[odd_val] + x[odd_val+1]
end
x = xsum
n = int(n/2)
end
return x[1]
end
# check its accuracy for a set logarithmically spaced n's. Since div2sum is slow,
# we won't go to very large n or use too many points
# restrict N to values of the form 2^n
n = 16
N = zeros(n)
for i = 1:n
N[i] = int(2^i)
end
err = Float64[]
err_perform = Float64[]
for n in N
x = rand(Float32, int(n))
xdouble = float64(x)
push!(err, abs(div2sum(x) - sum(xdouble)) / abs(sum(xdouble)))
push!(err_perform, abs(div2sum_performance(x) - sum(xdouble)) / abs(sum(xdouble)))
end
In [18]:
### make plots of divsum vs divsum_performance
using PyPlot
loglog(N, err, "bo-", label="div2sum")
loglog(N, err_perform, "ro--", label="div2sum_performance")
title("simple div2sum vs. performance")
xlabel("number of summands")
ylabel("relative error")
legend(loc="upper right",fancybox="true")
grid()
In [23]:
## compare speed
x = rand(Float32, 2^25)
@time div2sum(x)
@time div2sum_performance(x)
@time sum(x)
# show that they all give the same result
div2sum(x), div2sum_performance(x), sum(x)
Out[23]:
In [52]:
### show observed error vs. bound
n = 28
N = zeros(n)
for i = 1:n
N[i] = int(2^i)
end
err_perform = Float64[]
err_bound = zeros(length(N))
for i = 1:length(N)
n = N[i]
x = rand(Float32, int(n))
xdouble = float64(x)
## slope for error bound
err_bound[i] = (eps(Float32)*float32(Base.log2(n))*sum(abs(x)))
push!(err_perform, abs(sum(x) - sum(xdouble)) / abs(sum(xdouble)))
end
### make plots of divsum vs divsum_performance
using PyPlot
loglog(N, err_perform, "ro--",label="observed error")
loglog(N, err_bound, "bo-",label="error bound")
title("div2sum performance")
xlabel("number of summands")
ylabel("relative error")
legend(loc="upper right",fancybox="true")
grid()
In [ ]: