# functions for manipulating polynomials
# use this by pasting file contents into the code-area
# polynomials are assumed to be little-endian vectors of coefficients

# evaluate polynomial p at input x
pe(x, p) = sum[(i = 1), len(p), p[i] * x^(i - 1)]

# derivative of polynomial p (a new polynomial)
pd(p) = tail(map(p, (x, i => (i - 1) * x)))

# root of polynomial p near x (by Newton's method)
pr(p, x) = for[(dp = pd(p)), \
    10^-14 < abs(pe(x, p)), \
    (x = x - pe(x, p) / pe(x, dp)), \
    x]

# integral (antiderivative) of polynomial p (a new polynomial)
pa(p) = (.. 0) .. map(p, (x, i => /i * x))

# definite integral of polynomial p between bounds
pdi(p, a, b) = (q => pe(b, q) - pe(a, q))(pa(p))

# divide p by the trivial linear polynomial (-a, 1), ignoring remainder
dtl(p, a) = if[len(p) == 1, \
    zero(0), \
    (dtl((trim(trim(p)), p[len(p) - 1] + p[len(p)] * a), a), p[len(p)])]

# internal helper
discard(x) = 0

# all roots of polynomial p
pr(p) = if[len(p) < 2, \
    zero(0), \
    (r => (pr(dtl(p, r)), r))(pr(p, 0))]

# multiply polynomials p and q
pm(p, q) = sum[(k = 1), len(p), \
    (zero(k - 1) .. (q * p[k]) .. zero(len(p) - k))]

# inner product of polynomials over integration interval
ip(p, q, a, b) = pdi(pm(p, q), a, b)

# internal helper
chebygen(f) = n => pm((0, 2), f(n - 1)) - (f(n - 2), 0, 0)

# n-th Chebyshev polnomial of the first kind
chebycos(n) = chebygen(chebycos)(n)
chebycos(1) = (0, 1)
chebycos(0) = .. 1

# n-th Chebyshev polynomial of the second kind
chebysin(n) = chebygen(chebysin)(n)
chebysin(1) = (0, 2)
chebysin(0) = .. 1

# polynomial to an (integer) power
pp(p, n) = if[n / 2 == floor(n / 2), \
    (q => pm(q, q))(pp(p, n / 2)), \
    pm(p, pp(p, n - 1))]
pp(p, 0) = .. 1

# total of list of polynomials (of possible varying degree)
pt(l, n) = reduce(map(l, (p => p .. zero(n - len(p)))), (a, b => a + b))
pt(l) = pt(l, reduce(map(l, (p => len(p))), (a, b => if[a < b, b, a])))