Calc Lanczos polynomials

calc.py

# -*- coding: utf-8 -*-
#%%

from sympy import symbols, cos, pi, factor_list
from sympy.core.numbers import One
from sympy.utilities.lambdify import lambdify
from math import sin, pi
import numpy as np

def chebyshev(func, interval, degree=7, precision=20):
    n = degree + 1
    x, u = symbols('x u')
    a, b = interval
    x_to_u = (2 * x - a - b) / (b - a)
    u_to_x = (b - a) / 2 * u + (a + b) / 2
    chebyshev_nodes = cos((symbols('i') + 0.5) / n * pi)
    result_u = [ chebyshev_nodes.evalf(precision, subs={'i': i}) for i in range(n) ]
    result_x = [ u_to_x.evalf(precision, subs={u: i}) for i in result_u ]
    result_y = [ func(i) for i in result_x ]
    t = [One(), u]
    for _ in range(n-2):
        t.append(2 * u * t[-1] - t[-2])
    c = [ sum(result_y) / n ]
    for index in range(1, n):
        c.append( 2 * sum(t[index].evalf(precision, subs={u: i}) * j for i, j in zip(result_u, result_y)) / n )
    y = 1 * c[0]
    for i in range(1, n):
        y += t[i] * c[i]
    y = y.subs({u: x_to_u}).simplify()
    f = lambdify(x, y)
    f.formula = y
    return f
#%%

lanczos = lambda x: 2 * sin(pi*(x+0)) * sin(pi*(x+0)/2) / ((pi*(x+0))**2)

# That is the reversed list of polynomials
# Take rightmost polynomial to pop off the stack

def glsl_expr(varname, p):
    coeffs = p.as_poly().all_coeffs()
    expr = ""
    while len(coeffs) > 0:
        factor = coeffs[0]
        coeffs = coeffs[1:]
        if expr:
            expr = "(" + expr + ") * f + "
        expr += "{:0.4f}".format(factor)
    expr = expr.replace("+ -", "- ")
    print("vec4 {} = {};".format(varname, expr))

for degree in range(2, 5):
    print("/* Lanczos-{} */".format(degree))
    for interval in range(degree):
        for order in range(1, 6):
            lanczos = lambda x: degree * sin(pi*(x+interval)) * sin(pi*(x+interval)/degree) / ((pi*(x+interval))**2)
            f = chebyshev(lanczos, (0, 1), order)

            glsl_expr(
                "l{}_w{}_o{}".format(degree, interval, order),
                f.formula
            )

# %%

results.glsl

// All the kernels below take an input vec4 f which runs from 0 to 1 as sub-sample position

/* Lanczos-2 kernels for 0 > abs(x) > 1 */
vec4 l2_w0_o1 = (-1.1828) * f + 1.1298;
vec4 l2_w0_o2 = ((-0.2857) * f - 0.8025) * f + 1.0458;
vec4 l2_w0_o3 = (((1.5672) * f - 2.6445) * f + 0.0837) * f + 0.9976;
vec4 l2_w0_o4 = ((((-0.1262) * f + 1.7661) * f - 2.7215) * f + 0.0847) * f + 0.9983;
vec4 l2_w0_o5 = (((((-0.8556) * f + 2.0211) * f - 0.1210) * f - 2.0443) * f - 0.0003) * f + 1.0000;
/* Lanczos-2 kernels for 1 > abs(x) > 2 */
vec4 l2_w1_o1 = (0.0858) * f - 0.0792;
vec4 l2_w1_o2 = ((0.2379) * f - 0.1965) * f - 0.0249;
vec4 l2_w1_o3 = (((-0.7389) * f + 1.3652) * f - 0.6295) * f - 0.0004;
vec4 l2_w1_o4 = ((((0.3027) * f - 1.3178) * f + 1.7027) * f - 0.6892) * f + 0.0010;
vec4 l2_w1_o5 = (((((0.4259) * f - 0.7790) * f - 0.3527) * f + 1.3494) * f - 0.6436) * f + 0.0001;

/* Lanczos-3 */
vec4 l3_w0_o1 = (-1.1553) * f + 1.1305;
vec4 l3_w0_o2 = ((-0.4362) * f - 0.6392) * f + 1.0366;
vec4 l3_w0_o3 = (((1.3150) * f - 2.4044) * f + 0.0938) * f + 0.9972;
vec4 l3_w0_o4 = ((((0.0697) * f + 1.1386) * f - 2.2618) * f + 0.0556) * f + 0.9989;
vec4 l3_w0_o5 = (((((-0.5908) * f + 1.5478) * f - 0.1553) * f - 1.7997) * f - 0.0020) * f + 1.0000;
vec4 l3_w1_o1 = (0.0836) * f - 0.1080;
vec4 l3_w1_o2 = ((0.5461) * f - 0.5058) * f - 0.0187;
vec4 l3_w1_o3 = (((-0.7197) * f + 1.6203) * f - 0.9038) * f + 0.0023;
vec4 l3_w1_o4 = ((((-0.0846) * f - 0.5246) * f + 1.4755) * f - 0.8679) * f + 0.0008;
vec4 l3_w1_o5 = (((((0.4143) * f - 1.1236) * f + 0.3879) * f + 1.1483) * f - 0.8268) * f - 0.0000;
vec4 l3_w2_o1 = (-0.0287) * f + 0.0270;
vec4 l3_w2_o2 = ((-0.0955) * f + 0.0818) * f + 0.0073;
vec4 l3_w2_o3 = (((0.2538) * f - 0.4786) * f + 0.2263) * f - 0.0006;
vec4 l3_w2_o4 = ((((-0.0384) * f + 0.3183) * f - 0.5078) * f + 0.2287) * f - 0.0004;
vec4 l3_w2_o5 = (((((-0.1963) * f + 0.4574) * f - 0.1208) * f - 0.3487) * f + 0.2084) * f - 0.0000;

/* Lanczos-4 */
vec4 l4_w0_o1 = (-1.1448) * f + 1.1306;
vec4 l4_w0_o2 = ((-0.4894) * f - 0.5811) * f + 1.0333;
vec4 l4_w0_o3 = (((1.2200) * f - 2.3114) * f + 0.0952) * f + 0.9971;
vec4 l4_w0_o4 = ((((0.1273) * f + 0.9343) * f - 2.1056) * f + 0.0460) * f + 0.9991;
vec4 l4_w0_o5 = (((((-0.4983) * f + 1.3724) * f - 0.1541) * f - 1.7177) * f - 0.0023) * f + 1.0000;
vec4 l4_w1_o1 = (0.0781) * f - 0.1187;
vec4 l4_w1_o2 = ((0.6798) * f - 0.6421) * f - 0.0153;
vec4 l4_w1_o3 = (((-0.6703) * f + 1.6710) * f - 1.0035) * f + 0.0032;
vec4 l4_w1_o4 = ((((-0.2280) * f - 0.1914) * f + 1.3516) * f - 0.9337) * f + 0.0007;
vec4 l4_w1_o5 = (((((0.3670) * f - 1.1446) * f + 0.6094) * f + 1.0665) * f - 0.8982) * f - 0.0000;
vec4 l4_w2_o1 = (-0.0304) * f + 0.0433;
vec4 l4_w2_o2 = ((-0.2472) * f + 0.2328) * f + 0.0053;
vec4 l4_w2_o3 = (((0.2692) * f - 0.6435) * f + 0.3759) * f - 0.0018;
vec4 l4_w2_o4 = ((((0.1210) * f + 0.0139) * f - 0.4723) * f + 0.3381) * f - 0.0004;
vec4 l4_w2_o5 = (((((-0.2128) * f + 0.6535) * f - 0.4522) * f - 0.3059) * f + 0.3174) * f + 0.0000;
vec4 l4_w3_o1 = (0.0140) * f - 0.0133;
vec4 l4_w3_o2 = ((0.0507) * f - 0.0441) * f - 0.0033;
vec4 l4_w3_o3 = (((-0.1247) * f + 0.2379) * f - 0.1140) * f + 0.0004;
vec4 l4_w3_o4 = ((((0.0024) * f - 0.1231) * f + 0.2312) * f - 0.1109) * f + 0.0002;
vec4 l4_w3_o5 = (((((0.1035) * f - 0.2580) * f + 0.1066) * f + 0.1484) * f - 0.1004) * f + 0.0000;