LZX compressor improvements

[wimlib] / tools / log2_interpolation.r

diff --git a/tools/log2_interpolation.r b/tools/log2_interpolation.r
new file mode 100644 (file)
index 0000000..9523ec1
--- /dev/null
+++ b/tools/log2_interpolation.r
@@ -0,0 +1,47 @@
+#
+# This R language script finds a degree 2 interpolating polynomial for f =
+# log2(x) over the interval [1, 2).
+#
+
+f = function(x) { log2(x) }
+
+get.chebyshev.points = function(a, b, n)
+{
+    (a+b)/2 + (b-a)/2 * cos(((2*(1:n)-1)*pi)/(2*n))
+}
+
+build.vandermonde.matrix = function(x)
+{
+    n = length(x)
+    V = matrix(nrow=n, ncol=n)
+    for (j in 1:n)
+        V[,j] = x^(j-1)
+    return(V)
+}
+
+evaluate.polynomial = function(coeffs, x)
+{
+    y = coeffs[length(coeffs)]
+    for (i in (length(coeffs)-1):1)
+        y = y*x + coeffs[i]
+    return(y)
+}
+
+x.plot = seq(1, 2, length=1000)
+x.chebychev = get.chebyshev.points(1, 2, 3)
+V = build.vandermonde.matrix(x.chebychev)
+coeffs.a = solve(V, f(x.chebychev))
+coeffs.a = c(coeffs.a)
+polynomial.interp = function(x) { evaluate.polynomial(coeffs.a, x) }
+cat("Coefficients of degree 2 interpolating polynomial:\n")
+options(digits=10)
+cat(coeffs.a, "\n")
+pdf("polynomial-interp.pdf")
+plot(x.plot, f(x.plot), col="black", type="l", xlab="x", ylab="y",
+        main="f(x) and interpolating polynomial approximation")
+points(x.chebychev, f(x.chebychev), pch=19, col="red")
+lines(x.plot, polynomial.interp(x.plot), col="blue")
+legend("topleft", pch=c(NA, NA, 19),
+          col=c("black", "blue", "red"), lty=c(1,1,0),
+       legend=c("f(x)", "Interpolating polynomial", "Chebychev points"),
+       inset=0.07)