LZX compressor improvements
[wimlib] / tools / log2_interpolation.r
1 #
2 # This R language script finds a degree 2 interpolating polynomial for f =
3 # log2(x) over the interval [1, 2).
4 #
5
6 f = function(x) { log2(x) }
7
8 get.chebyshev.points = function(a, b, n)
9 {
10     (a+b)/2 + (b-a)/2 * cos(((2*(1:n)-1)*pi)/(2*n))
11 }
12
13 build.vandermonde.matrix = function(x)
14 {
15     n = length(x)
16     V = matrix(nrow=n, ncol=n)
17     for (j in 1:n)
18         V[,j] = x^(j-1)
19     return(V)
20 }
21
22 evaluate.polynomial = function(coeffs, x)
23 {
24     y = coeffs[length(coeffs)]
25     for (i in (length(coeffs)-1):1)
26         y = y*x + coeffs[i]
27     return(y)
28 }
29
30 x.plot = seq(1, 2, length=1000)
31 x.chebychev = get.chebyshev.points(1, 2, 3)
32 V = build.vandermonde.matrix(x.chebychev)
33 coeffs.a = solve(V, f(x.chebychev))
34 coeffs.a = c(coeffs.a)
35 polynomial.interp = function(x) { evaluate.polynomial(coeffs.a, x) }
36 cat("Coefficients of degree 2 interpolating polynomial:\n")
37 options(digits=10)
38 cat(coeffs.a, "\n")
39 pdf("polynomial-interp.pdf")
40 plot(x.plot, f(x.plot), col="black", type="l", xlab="x", ylab="y",
41          main="f(x) and interpolating polynomial approximation")
42 points(x.chebychev, f(x.chebychev), pch=19, col="red")
43 lines(x.plot, polynomial.interp(x.plot), col="blue")
44 legend("topleft", pch=c(NA, NA, 19),
45            col=c("black", "blue", "red"), lty=c(1,1,0),
46        legend=c("f(x)", "Interpolating polynomial", "Chebychev points"),
47        inset=0.07)