Skip to content

Standard

Polynomials in standard monomial basis.

\[ p(x)=\sum_{k=0}^na_kx^k \qquad P=(x_n)_{n\in\mathbb{N}_0}=(1, x, x^2, \dots) \]

Prefixed by poly... (polynomial).

All functions accept single exhaustible iterables, if not stated otherwise.

If the result of a function is a polynomial, it is returned as a tuple.

Associative operations, like polyadd and polymul, allow arbitrary many arguments, even none.

The functions are type-independent. However, the data types used must support necessary scalar operations. For instance, for polynomial addition, components must be addable.

Docstring convention

Summary

Math notation

Complexity

For a polynomial of degree \(n\) there will be

  • \(x\) scalar additions (add),
  • \(y\) scalar subtractions, ...

Order: pos, neg, add, sub, mul

Notes

design choices

See also

other implementations, delegations/wrappings, ...

References

numpy equivalents, Wikipedia, ...)

creation

polyzero = ()

Zero polynomial.

\[ 0 \qquad 0^{(P)}=\vec{0} \]

An empty tuple: ().

Notes

Why give the zero polynomial a distinguished representation (not just [0] like numpy.polynomial)?

The additive neutral element seems worth handling exceptionally. It is mathematically different (different degree) and results in functions like polymul being more time and memory efficient.

See also
References

polyone = (1,)

Constant one polynomial.

\[ 1 \qquad 1^{(P)}=\begin{pmatrix} 1 \end{pmatrix} \]

A tuple with a single one: (1,).

See also
References

polyx = (0, 1)

Identity polynomial.

\[ x \qquad x^{(P)}=\begin{pmatrix} 0 \\ 1 \end{pmatrix} \]

A tuple with a zero and a one: (0, 1).

See also
References

polymono(n, c=1, zero=0)

Return a monomial of degree n.

\[ cx^n \qquad (cx^n)^{(P)}=c\vec{e}_n \]

Returns a tuple with n zeros followed by c or polyzero if \(n<0\).

See also
Source code in poly\standard\creation.py 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def polymono(n, c=1, zero=0):
    r"""Return a monomial of degree `n`.

    $$
        cx^n \qquad (cx^n)^{(P)}=c\vec{e}_n
    $$

    Returns a tuple with `n` zeros followed by `c` or [`polyzero`][poly.standard.polyzero] if $n<0$.

    See also
    --------
    - constants: [`polyzero`][poly.standard.polyzero], [`polyone`][poly.standard.polyone], [`polyx`][poly.standard.polyx]
    - for all monomials: [`polymonos`][poly.standard.polymonos]
    - wraps: [`vector.vecbasis`](https://goessl.github.io/vector/functional/#vector.functional.creation.vecbasis)
    """
    return vecbasis(n, c=c, zero=zero) if n>=0 else polyzero

polymonos(start=0, c=1, zero=0)

Yield all monomials.

\[ \left(x^n\right)_\mathbb{n\in\mathbb{N_0}} = \left(1, x, x^2, \dots \right) \]
See also
Source code in poly\standard\creation.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def polymonos(start=0, c=1, zero=0):
    r"""Yield all monomials.

    $$
        \left(x^n\right)_\mathbb{n\in\mathbb{N_0}} = \left(1, x, x^2, \dots \right)
    $$

    See also
    --------
    - for single monomial: [`polymono`][poly.standard.polymono]
    - wraps: [`vector.vecbases`](https://goessl.github.io/vector/functional/#vector.functional.creation.vecbases)
    """
    yield from vecbases(start=start, c=c, zero=zero)

polyrand(n)

Return a random polynomial of degree n.

\[ \sum_{k=0}^na_kx^k \quad a_k\sim\mathcal{U}([0, 1[) \]

The coefficients are sampled from a uniform distribution in [0, 1[.

See also
Source code in poly\standard\creation.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def polyrand(n):
    r"""Return a random polynomial of degree `n`.

    $$
        \sum_{k=0}^na_kx^k \quad a_k\sim\mathcal{U}([0, 1[)
    $$

    The coefficients are sampled from a uniform distribution in `[0, 1[`.

    See also
    --------
    - wraps: [`vector.vecrand`](https://goessl.github.io/vector/functional/#vector.functional.creation.vecrand)
    """
    return vecrand(n+1)

polyrandn(n, normed=True, mu=0, sigma=1)

Return a random polynomial of degree n.

\[ \sum_{k=0}^na_kx^k \quad a_k\sim\mathcal{N}(\mu, \sigma) \]

The coefficients are sampled from a normal distribution.

Normed with respect to the euclidian vector norm \(\sum_k|a_k|^2\).

See also
Source code in poly\standard\creation.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def polyrandn(n, normed=True, mu=0, sigma=1):
    r"""Return a random polynomial of degree `n`.

    $$
        \sum_{k=0}^na_kx^k \quad a_k\sim\mathcal{N}(\mu, \sigma)
    $$

    The coefficients are sampled from a normal distribution.

    Normed with respect to the euclidian vector norm $\sum_k|a_k|^2$.

    See also
    --------
    - wraps: [`vector.vecrandn`](https://goessl.github.io/vector/functional/#vector.functional.creation.vecrandn)
    """
    return vecrandn(n+1, normed=normed, mu=mu, sigma=sigma)

polyfromroots(*xs, one=1)

Return the polynomial with the given roots.

\[ \prod_k(x-x_k) \]
Complexity

For \(n\) roots there will be

  • \(n\) scalar negations (neg),
  • \(\frac{n(n-1)}{2}\) scalar additions (add) &
  • \(\begin{cases}(n+2)(n-1)&n\ge1\\0&n\le1\end{cases}\) scalar multiplications (mul).
References
Source code in poly\standard\creation.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def polyfromroots(*xs, one=1):
    r"""Return the polynomial with the given roots.

    $$
        \prod_k(x-x_k)
    $$

    Complexity
    ----------
    For $n$ roots there will be

    - $n$ scalar negations (`neg`),
    - $\frac{n(n-1)}{2}$ scalar additions (`add`) &
    - $\begin{cases}(n+2)(n-1)&n\ge1\\0&n\le1\end{cases}$ scalar multiplications (`mul`).

    References
    ----------
    - Recipe: [more_itertools.polynomial_from_roots](https://more-itertools.readthedocs.io/en/stable/api.html#more_itertools.polynomial_from_roots)
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polyfromroots`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyfromroots.html)
    """
    #return polymul(*zip(map(neg, xs), repeat(one)), method='naive', one=one)
    r = (one, )
    for x in xs:
        r = polysub(polymulx(r), polyscalarmul(x, r))
    return r

utility

polydeg(p)

Return the degree of a polynomial.

\[ \deg(p) \]

Doesn't handle leading zeros, use polytrim if needed.

\(\deg(0)=-1\) is used for the empty zero polynomial.

Notes

\(\deg(0)=-\infty\) is more commonly used but the expected return type is an int where -math.inf is of type float. Therefore the \(\deg(0)=-1\) convention was choosen to keep the return type consistent.

See also
References
Source code in poly\standard\utility.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def polydeg(p):
    r"""Return the degree of a polynomial.

    $$
        \deg(p)
    $$

    Doesn't handle leading zeros, use [`polytrim`][poly.standard.polytrim]
    if needed.

    $\deg(0)=-1$ is used for the empty zero polynomial.

    Notes
    -----
    $\deg(0)=-\infty$ is more commonly used but the expected return type is an
    `int` where `-math.inf` is of type `float`. Therefore the $\deg(0)=-1$
    convention was choosen to keep the return type consistent.

    See also
    --------
    - wraps: [`vector.veclen`](https://goessl.github.io/vector/functional/#vector.functional.utility.veclen)

    References
    ----------
    - [Wikipedia - Degree of a polynomial - Degree of the zero polynomial](https://en.wikipedia.org/wiki/Degree_of_a_polynomial#Degree_of_the_zero_polynomial)
    """
    return veclen(p) - 1

polyeq(p, q)

Return if two polynomials are equal.

\[ p \overset{?}{=} q \]
Complexity

For two polynomials of degrees \(n\) & \(m\) there will be at most

  • \(\min\{n, m\}+1\) scalar comparisons (eq) &
  • \(|n-m|\) scalar boolean evaluations (bool).
See also
Source code in poly\standard\utility.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def polyeq(p, q):
    r"""Return if two polynomials are equal.

    $$
        p \overset{?}{=} q
    $$

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be at most

    - $\min\{n, m\}+1$ scalar comparisons (`eq`) &
    - $|n-m|$ scalar boolean evaluations (`bool`).

    See also
    --------
    - wraps: [`vector.veceq`](https://goessl.github.io/vector/functional/#vector.functional.utility.veceq)
    """
    return veceq(p, q)

polytrim(p, tol=1e-09)

Remove all leading near zero (abs(a_i)<=tol) coefficients.

\[ \sum_{k=0}^na_kx^k \ \text{where} \ n=\max\{\, k\mid |a_k|>\text{tol}\,\}\cup\{-1\} \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar absolute evaluations (abs) &
  • \(n+1\) scalar comparisons (gt).
See also
Source code in poly\standard\utility.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def polytrim(p, tol=1e-9):
    r"""Remove all leading near zero (`abs(a_i)<=tol`) coefficients.

    $$
        \sum_{k=0}^na_kx^k \ \text{where} \ n=\max\{\, k\mid |a_k|>\text{tol}\,\}\cup\{-1\}
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar absolute evaluations (`abs`) &
    - $n+1$ scalar comparisons (`gt`).

    See also
    --------
    - wraps: [`vector.vectrim`](https://goessl.github.io/vector/functional/#vector.functional.utility.vectrim)
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polytrim`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polytrim.html)
    """
    return vectrim(p, tol=tol)

evaluation

polyval(p, x, method='horner')

Return the value of polynomial p evaluated at point x.

\[ p(x) \]

Available methods are

See also
References
Source code in poly\standard\evaluation.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def polyval(p, x, method='horner'):
    """Return the value of polynomial `p` evaluated at point `x`.

    $$
        p(x)
    $$

    Available methods are

    - [`naive`][poly.standard.polyval_naive],
    - [`iterative`][poly.standard.polyval_iterative] &
    - [`horner`][poly.standard.polyval_horner] (`p` must be reversible).

    See also
    --------
    - implementations: [`polyval_naive`][poly.standard.polyval_naive],
    [`polyval_iterative`][poly.standard.polyval_iterative],
    [`polyval_horner`][poly.standard.polyval_horner]
    - for consecutive monomials: [`polyvals`][poly.standard.polyvals]
    - for $x=0$: [`polyvalzero`][poly.standard.polyvalzero]
    - for polynomial arguments: [`polycom`][poly.standard.polycom]

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polyval`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyval.html)
    """
    match method:
        case 'naive':
            return polyval_naive(p, x)
        case 'iterative':
            return polyval_iterative(p, x)
        case 'horner':
            return polyval_horner(p, x)
        case _:
            raise ValueError('Invalid method')

polyval_naive(p, x)

Return the value of polynomial p evaluated at point x.

\[ p(x) \]

Uses naive repeated multiplication to calculate monomials individually.

See also
References
Source code in poly\standard\evaluation.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
def polyval_naive(p, x):
    """Return the value of polynomial `p` evaluated at point `x`.

    $$
        p(x)
    $$

    Uses naive repeated multiplication to calculate monomials individually.

    See also
    --------
    - for any implementation: [`polyval`][poly.standard.polyval]
    - other implementations:
    [`polyval_iterative`][poly.standard.polyval_iterative],
    [`polyval_horner`][poly.standard.polyval_horner]
    - for polynomial arguments: [`polycom_naive`][poly.standard.polycom_naive]

    References
    ----------
    - [more_itertools.polynomial_eval](https://more-itertools.readthedocs.io/en/stable/api.html#more_itertools.polynomial_eval)
    - [Wikipedia - Polynomial evaluation](https://en.wikipedia.org/wiki/Polynomial_evaluation)
    - [Wikipedia - Horner's method - Efficiency](https://en.wikipedia.org/wiki/Horner%27s_method#Efficiency)
    """
    p = iter(p)
    a0 = next(p, type(x)(0))
    return sumprod_default(p, chain((x,), map(pow, repeat(x), count(start=2))), initial=a0, default=MISSING)

polyval_iterative(p, x)

Return the value of polynomial p evaluated at point x.

\[ p(x) \]

Uses iterative multiplication to calculate monomials consecutively.

Complexity

For a polynomial of degree \(n\) there will be

  • \(\begin{cases}n&n\ge0\\0&n\le0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}2n-1&n>0\\0&n\le0\end{cases}\) scalar multiplications (mul).
See also
References
Source code in poly\standard\evaluation.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def polyval_iterative(p, x):
    r"""Return the value of polynomial `p` evaluated at point `x`.

    $$
        p(x)
    $$

    Uses iterative multiplication to calculate monomials consecutively.

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $\begin{cases}n&n\ge0\\0&n\le0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}2n-1&n>0\\0&n\le0\end{cases}$ scalar multiplications (`mul`).

    See also
    --------
    - for any implementation: [`polyval`][poly.standard.polyval]
    - other implementations:
    [`polyval_naive`][poly.standard.polyval_naive],
    [`polyval_horner`][poly.standard.polyval_horner]
    - uses: [`polyvals`][poly.standard.polyvals]
    - for polynomial arguments: [`polycom_iterative`][poly.standard.polycom_iterative]

    References
    ----------
    - [Wikipedia - Horner's method - Efficiency](https://en.wikipedia.org/wiki/Horner%27s_method#Efficiency)
    """
    p = iter(p)
    a0 = next(p, type(x)(0))
    return sumprod_default(p, polyvals(x, start=1), initial=a0, default=MISSING)

polyval_horner(p, x)

Return the value of polynomial p evaluated at point x.

\[ p(x) \]

Uses Horner's method.

p must be reversible.

Complexity

For a polynomial of degree \(n\) there will be

  • \(\begin{cases}n&n\ge0\\0&n\le0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}n&n\ge0\\0&n\le0\end{cases}\) scalar multiplications (mul).
See also
References
Source code in poly\standard\evaluation.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
def polyval_horner(p, x):
    r"""Return the value of polynomial `p` evaluated at point `x`.

    $$
        p(x)
    $$

    Uses Horner's method.

    `p` must be reversible.

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $\begin{cases}n&n\ge0\\0&n\le0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}n&n\ge0\\0&n\le0\end{cases}$ scalar multiplications (`mul`).

    See also
    --------
    - for any implementation: [`polyval`][poly.standard.polyval]
    - other implementations:
    [`polyval_naive`][poly.standard.polyval_naive],
    [`polyval_iterative`][poly.standard.polyval_iterative]
    - for polynomial arguments: [`polycom_horner`][poly.standard.polycom_horner]

    References
    ----------
    - [Wikipedia - Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method)
    """
    p = iter(reversed(p))
    an = next(p, type(x)(0))
    return reduce_default(lambda a, pi: a*x+pi, p, initial=an, default=MISSING)

polyvals(x, start=0)

Yield the powers of the value x.

\[ (x^n)_{n\in\mathbb{N}_0} = (1, x, x^2, \dots) \]

Uses iterative multiplication to calculate monomials consecutively.

See also
Source code in poly\standard\evaluation.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def polyvals(x, start=0):
    r"""Yield the powers of the value `x`.

    $$
        (x^n)_{n\in\mathbb{N}_0} = (1, x, x^2, \dots)
    $$

    Uses iterative multiplication to calculate monomials consecutively.

    See also
    --------
    - used by: [`polyval_iterative`][poly.standard.polyval_iterative]
    - for polynomial arguments: [`polypows`][poly.standard.polypows]
    """
    if start <= 0:
        yield type(x)(1)
    yield (y := prod_default(repeat(x, start), initial=MISSING, default=x))
    while True:
        y *= x
        yield y

polyvalzero(p, zero=0)

Return the value of polynomial p evaluated at point 0.

\[ p(0) \]

More efficient than polyval(p, 0).

Complexity

There are no scalar arithmetic operations.

See also
Source code in poly\standard\evaluation.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def polyvalzero(p, zero=0):
    """Return the value of polynomial `p` evaluated at point 0.

    $$
        p(0)
    $$

    More efficient than `polyval(p, 0)`.

    Complexity
    ----------
    There are no scalar arithmetic operations.

    See also
    --------
    - for any argument: [`polyval`][poly.standard.polyval]
    """
    return next(iter(p), zero)

polycom(p, q, method='iterative')

Return the polynomial composition of p & q.

\[ p\circ q \]

q must be a sequence.

Available methods are

See also
Source code in poly\standard\evaluation.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def polycom(p, q, method='iterative'):
    r"""Return the polynomial composition of `p` & `q`.

    $$
        p\circ q
    $$

    `q` must be a sequence.

    Available methods are

    - [`naive`][poly.standard.polycom_naive],
    - [`iterative`][poly.standard.polycom_iterative] &
    - [`horner`][poly.standard.polycom_horner] (`p` must be reversible).

    See also
    --------
    - implementations: [`polycom_naive`][poly.standard.polycom_naive],
    [`polycom_iterative`][poly.standard.polycom_iterative],
    [`polycom_horner`][poly.standard.polycom_horner]
    - for $q=x-s$: [`polyshift`][poly.standard.polyshift]
    - for scalar arguments: [`polyval`][poly.standard.polyval]
    """
    match method:
        case 'naive':
            return polycom_naive(p, q)
        case 'iterative':
            return polycom_iterative(p, q)
        case 'horner':
            return polycom_horner(p, q)
        case _:
            raise ValueError('Invalid method')

polycom_naive(p, q)

Return the polynomial composition of p & q.

\[ p\circ q \]

Uses naive repeated multiplication to calculate monomials individually.

q must be a sequence.

Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\begin{cases}\frac{n^3m^2}{6}+\frac{n^2m}{2}-\frac{nm^2}{6}-\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}\frac{n^3m^2}{6}+\frac{n^3m}{6}+n^2m-\frac{nm^2}{6}-\frac{nm}{6}+\frac{n^2}{2}+\frac{n}{2}&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}\) scalar multiplications (´mul´).
See also
Source code in poly\standard\evaluation.py 
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def polycom_naive(p, q):
    r"""Return the polynomial composition of `p` & `q`.

    $$
        p\circ q
    $$

    Uses naive repeated multiplication to calculate monomials individually.

    `q` must be a sequence.

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\begin{cases}\frac{n^3m^2}{6}+\frac{n^2m}{2}-\frac{nm^2}{6}-\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}\frac{n^3m^2}{6}+\frac{n^3m}{6}+n^2m-\frac{nm^2}{6}-\frac{nm}{6}+\frac{n^2}{2}+\frac{n}{2}&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}$ scalar multiplications (´mul´).

    See also
    --------
    - for any implementation: [`polycom`][poly.standard.polycom]
    - other implementations:
    [`polycom_iterative`][poly.standard.polycom_iterative],
    [`polycom_horner`][poly.standard.polycom_horner]
    - for scalar arguments: [`polyval_naive`][poly.standard.polyval_naive]
    """
    p = iter(p)
    a0 = tuple(islice(p, 1)) #(a0,) or polyzero
    return polyadd(a0, *(polyscalarmul(pi, polypow_naive(q, i)) for i, pi in enumerate(p, start=1)))

polycom_iterative(p, q)

Return the polynomial composition of p & q.

\[ p\circ q \]

Uses iterative multiplication to calculate monomials consecutively.

q must be a sequence.

Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\begin{cases}\frac{n^2m^2}{2}+\frac{n^2m}{2}-\frac{nm^2}{2}-\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}\frac{n^2m^2}{2}+n^2m-\frac{nm^2}{2}+nm+2n-m-1&n\ge1\land m\ge0\\0&n\le1\lor m\le0\end{cases}\) scalar multiplications (´mul´).
See also
Source code in poly\standard\evaluation.py 
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def polycom_iterative(p, q):
    r"""Return the polynomial composition of `p` & `q`.

    $$
        p\circ q
    $$

    Uses iterative multiplication to calculate monomials consecutively.

    `q` must be a sequence.

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\begin{cases}\frac{n^2m^2}{2}+\frac{n^2m}{2}-\frac{nm^2}{2}-\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}\frac{n^2m^2}{2}+n^2m-\frac{nm^2}{2}+nm+2n-m-1&n\ge1\land m\ge0\\0&n\le1\lor m\le0\end{cases}$ scalar multiplications (´mul´).

    See also
    --------
    - for any implementation: [`polycom`][poly.standard.polycom]
    - other implementations:
    [`polycom_naive`][poly.standard.polycom_naive],
    [`polycom_horner`][poly.standard.polycom_horner]
    - uses: [`polypows`][poly.standard.arithmetic.polypows]
    - for scalar arguments: [`polyval_iterative`][poly.standard.polyval_iterative]
    """
    p = iter(p)
    a0 = tuple(islice(p, 1))
    return polyadd(a0, *map(polyscalarmul, p, polypows(q, start=1)))

polycom_horner(p, q)

Return the polynomial composition of p & q.

\[ p\circ q \]

Uses Horner's method.

q must be reversible.

Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\begin{cases}0&m\ge0\lor n<0\\n&m<0\land n>=0\end{cases}\) scalar unary pluses (pos),
  • \(\begin{cases}\frac{n^2m^2}{2}-\frac{nm^2}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}\frac{n^2m^2}{2}+\frac{n^2m}{2}-\frac{nm^2}{2}+\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}\) scalar multiplications (´mul´).
See also
Source code in poly\standard\evaluation.py 
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
def polycom_horner(p, q):
    r"""Return the polynomial composition of `p` & `q`.

    $$
        p\circ q
    $$

    Uses Horner's method.

    `q` must be reversible.

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\begin{cases}0&m\ge0\lor n<0\\n&m<0\land n>=0\end{cases}$ scalar unary pluses (`pos`),
    - $\begin{cases}\frac{n^2m^2}{2}-\frac{nm^2}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}\frac{n^2m^2}{2}+\frac{n^2m}{2}-\frac{nm^2}{2}+\frac{nm}{2}+n&n\ge0\land m\ge0\\0&n\le0\lor m\le0\end{cases}$ scalar multiplications (´mul´).

    See also
    --------
    - for any implementation: [`polycom`][poly.standard.polycom]
    - other implementations:
    [`polycom_naive`][poly.standard.polycom_naive],
    [`polycom_iterative`][poly.standard.polycom_iterative]
    - for scalar arguments: [`polyval_horner`][poly.standard.polyval_horner]
    """
    p = iter(reversed(p))
    an = tuple(islice(p, 1))
    return reduce_default(lambda a, pi: polyaddc(polymul_naive(a, q), pi), p, initial=an, default=MISSING)

polyshift(p, s, one=1)

Return the polynomial p shifted by s on the abscissa.

\[ p(x - s) \]

TODO: https://math.stackexchange.com/a/694571/1170417

See also
Source code in poly\standard\evaluation.py 
309
310
311
312
313
314
315
316
317
318
319
320
321
322
def polyshift(p, s, one=1):
    """Return the polynomial `p` shifted by `s` on the abscissa.

    $$
        p(x - s)
    $$

    TODO: https://math.stackexchange.com/a/694571/1170417

    See also
    --------
    - for polynomial argument: [`polycom`][poly.standard.polycom]
    """
    return polycom_naive(p, (-s, one))

polyscale(p, a)

Return the polynomial p scaled by 1/a on the abscissa.

\[ p(ax) \]

More efficient than polycom(p, polyscalarmul(a, polyx)).

See also
Source code in poly\standard\evaluation.py 
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
def polyscale(p, a):
    """Return the polynomial `p` scaled by 1/`a` on the abscissa.

    $$
        p(ax)
    $$

    More efficient than `polycom(p, polyscalarmul(a, polyx))`.

    See also
    --------
    - for polynomial argument: [`polycom`][poly.standard.polycom]
    """
    p = iter(p)
    return tuple(chain(islice(p, 1), veclhadamard(p, polyvals(a, start=1))))

arithmetic

polypos(p)

Return the polynomial with the unary positive operator applied.

\[ +p \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar unary plus operations (pos).
See also
Source code in poly\standard\arithmetic.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def polypos(p):
    """Return the polynomial with the unary positive operator applied.

    $$
        +p
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar unary plus operations (`pos`).

    See also
    --------
    - wraps: [`vector.vecpos`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecpos)
    """
    return vecpos(p)

polyneg(p)

Return the polynomial with the unary negative operator applied.

\[ -p \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar negations (neg).
See also
Source code in poly\standard\arithmetic.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def polyneg(p):
    """Return the polynomial with the unary negative operator applied.

    $$
        -p
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar negations (`neg`).

    See also
    --------
    - wraps: [`vector.vecneg`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecneg)
    """
    return vecneg(p)

polyadd(*ps)

Return the sum of polynomials.

\[ p_0 + p_1 + \cdots \]
Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\min\{n, m\}+1\) scalar additions (add).
See also
References
Source code in poly\standard\arithmetic.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def polyadd(*ps):
    r"""Return the sum of polynomials.

    $$
        p_0 + p_1 + \cdots
    $$

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\min\{n, m\}+1$ scalar additions (`add`).

    See also
    --------
    - for constant or monomial summand: [`polyaddc`][poly.standard.polyaddc]
    - wraps: [`vector.vecadd`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecadd)

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polyadd`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyadd.html)
    """
    return vecadd(*ps)

polyaddc(p, c, n=0)

Return the sum of polynomial p and a monomial of degree n.

\[ p + cx^n \]

More efficient than polyadd(p, polymono(n, c)).

Complexity

There will be

  • one scalar addition (add) if \(n\le\deg p\) or
  • one unary plus operations (pos) otherwise.
See also
Source code in poly\standard\arithmetic.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def polyaddc(p, c, n=0):
    r"""Return the sum of polynomial `p` and a monomial of degree `n`.

    $$
        p + cx^n
    $$

    More efficient than `polyadd(p, polymono(n, c))`.

    Complexity
    ----------
    There will be

    - one scalar addition (`add`) if $n\le\deg p$ or
    - one unary plus operations (`pos`) otherwise.

    See also
    --------
    - for polynomial summand: [`polyadd`][poly.standard.polyadd]
    - wraps: [`vector.vecaddc`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecaddc)
    """
    return vecaddc(p, c, i=n)

polysub(p, q)

Return the difference of two polynomials.

\[ p - q \]
Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\min\{n, m\}+1\) scalar subtractions (sub) &
  • \(\begin{cases}m-n&m\ge n\\0&m\le n\end{cases}\) negations (neg).
See also
References
Source code in poly\standard\arithmetic.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def polysub(p, q):
    r"""Return the difference of two polynomials.

    $$
        p - q
    $$

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\min\{n, m\}+1$ scalar subtractions (`sub`) &
    - $\begin{cases}m-n&m\ge n\\0&m\le n\end{cases}$ negations (`neg`).

    See also
    --------
    - for constant or monomial subtrahend: [`polysubc`][poly.standard.polysubc]
    - wraps: [`vector.vecsub`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecsub)

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polysub`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polysub.html)
    """
    return vecsub(p, q)

polysubc(p, c, n=0)

Return the difference of polynomial p and a monomial of degree n.

\[ p - cx^n \]

More efficient than polysub(p, polymono(n, c)).

Complexity

There will be

  • one scalar subtraction (sub) if \(n\le\deg p\) or
  • one scalar negation (neg) otherwise.
See also
Source code in poly\standard\arithmetic.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def polysubc(p, c, n=0):
    r"""Return the difference of polynomial `p` and a monomial of degree `n`.

    $$
        p - cx^n
    $$

    More efficient than `polysub(p, polymono(n, c))`.

    Complexity
    ----------
    There will be

    - one scalar subtraction (`sub`) if $n\le\deg p$ or
    - one scalar negation (`neg`) otherwise.

    See also
    --------
    - for polynomial subtrahend: [`polysub`][poly.standard.polysub]
    - wraps: [`vector.vecsubc`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecsubc)
    """
    return vecsubc(p, c, i=n)

polyscalarmul(a, p)

Return the product of a scalar and a polynomial.

\[ a\cdot p \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar multiplications (rmul).
See also
Source code in poly\standard\arithmetic.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def polyscalarmul(a, p):
    r"""Return the product of a scalar and a polynomial.

    $$
        a\cdot p
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar multiplications (`rmul`).

    See also
    --------
    - for polynomial factor: [`polymul`][poly.standard.polymul]
    - wraps: [`vector.vecmul`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecmul)
    """
    return vecmul(a, p)

polyscalartruediv(p, a)

Return the true division of a polynomial and a scalar.

\[ \frac{p}{a} \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar true divisions (truediv).
See also
Source code in poly\standard\arithmetic.py 
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
def polyscalartruediv(p, a):
    r"""Return the true division of a polynomial and a scalar.

    $$
        \frac{p}{a}
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar true divisions (`truediv`).

    See also
    --------
    - wraps: [`vector.vectruediv`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vectruediv)
    """
    return vectruediv(p, a)

polyscalarfloordiv(p, a)

Return the floor division of a polynomial and a scalar.

\[ \sum_k\left\lfloor\frac{a_k}{a}\right\rfloor x^k \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar floor divisions (floordiv).
See also
Source code in poly\standard\arithmetic.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def polyscalarfloordiv(p, a):
    r"""Return the floor division of a polynomial and a scalar.

    $$
        \sum_k\left\lfloor\frac{a_k}{a}\right\rfloor x^k
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar floor divisions (`floordiv`).

    See also
    --------
    - included in: [`polyscalardivmod`][poly.standard.polyscalardivmod]
    - wraps: [`vector.vecfloordiv`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecfloordiv)
    """
    return vecfloordiv(p, a)

polyscalarmod(p, a)

Return the elementwise mod of a polynomial and a scalar.

\[ \sum_k\left(a_k \mod a \right)x^k \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar modulos (mod).
See also
Source code in poly\standard\arithmetic.py 
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
def polyscalarmod(p, a):
    r"""Return the elementwise mod of a polynomial and a scalar.

    $$
        \sum_k\left(a_k \mod a \right)x^k
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar modulos (`mod`).

    See also
    --------
    - included in: [`polyscalardivmod`][poly.standard.polyscalardivmod]
    - wraps: [`vector.vecmod`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecmod)
    """
    return vecmod(p, a)

polyscalardivmod(p, a)

Return the elementwise divmod of a polynomial and a scalar.

\[ \sum_k\left\lfloor\frac{a_k}{a}\right\rfloor x^k, \ \sum_k\left(a_k \mod a \right)x^k \]
Complexity

For a polynomial of degree \(n\) there will be

  • \(n+1\) scalar divmods (divmod).
See also
Source code in poly\standard\arithmetic.py 
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
def polyscalardivmod(p, a):
    r"""Return the elementwise divmod of a polynomial and a scalar.

    $$
        \sum_k\left\lfloor\frac{a_k}{a}\right\rfloor x^k, \ \sum_k\left(a_k \mod a \right)x^k
    $$

    Complexity
    ----------
    For a polynomial of degree $n$ there will be

    - $n+1$ scalar divmods (`divmod`).

    See also
    --------
    - combines: [`polyscalarfloordiv`][poly.standard.polyscalarfloordiv] & [`polyscalarmod`][poly.standard.polyscalarmod]
    - wraps: [`vector.vecdivmod`](https://goessl.github.io/vector/functional/#vector.functional.vector_space.vecdivmod)
    """
    return vecdivmod(p, a)

polymul(*ps, method='naive', one=1)

Return the product of polynomials.

\[ p_0 p_1 \cdots \]

Available methods are

See also
References
Source code in poly\standard\arithmetic.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
def polymul(*ps, method='naive', one=1):
    r"""Return the product of polynomials.

    $$
        p_0 p_1 \cdots
    $$

    Available methods are

    - [`naive`][poly.standard.polymul_naive] &
    - [`karatsuba`][poly.standard.polymul_karatsuba].

    See also
    --------
    - implementations: [`polymul_naive`][poly.standard.arithmetic.polymul_naive],
    [`polymul_karatsuba`][poly.standard.arithmetic.polymul_karatsuba]
    - for scalar factor: [`polyscalarmul`][poly.standard.arithmetic.polyscalarmul]
    - for monomial factor: [`polymulx`][poly.standard.arithmetic.polymulx]

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polymul`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polymul.html)
    """
    ps = tuple(map(tuple, ps))
    match method:
        case 'naive':
            return reduce_default(polymul_naive, ps, default=(one,))
        case 'karatsuba':
            return reduce_default(polymul_karatsuba, ps, default=(one,))
        case _:
            raise ValueError('Invalid method')

polymul_naive(p, q)

Return the product of two polynomials.

\[ pq \]

Uses naive multiplication and summation.

q must be a sequence.

Complexity

For two polynomials of degrees \(n\) & \(m\) there will be

  • \(\begin{cases}nm&n\ge1\land m\ge1\\0&n\le0\lor m\le0\end{cases}\) scalar additions (add) &
  • \((n+1)(m+1)\) scalar mutliplications (mul).
See also
Source code in poly\standard\arithmetic.py 
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
def polymul_naive(p, q):
    r"""Return the product of two polynomials.

    $$
        pq
    $$

    Uses naive multiplication and summation.

    `q` must be a sequence.

    Complexity
    ----------
    For two polynomials of degrees $n$ & $m$ there will be

    - $\begin{cases}nm&n\ge1\land m\ge1\\0&n\le0\lor m\le0\end{cases}$ scalar additions (`add`) &
    - $(n+1)(m+1)$ scalar mutliplications (`mul`).

    See also
    --------
    - for any implementation: [`polymul`][poly.standard.arithmetic.polymul]
    - other implementations:
    [`polymul_karatsuba`][poly.standard.arithmetic.polymul_karatsuba]
    """
    if not q:
        return () #polyzero

    r, sentinel = [], object()
    for i, pi in enumerate(p):
        r.extend([sentinel] * (i+len(q) - len(r)))
        for j, qj in enumerate(q):
            if r[i+j] is sentinel:
                r[i+j] = pi * qj
            else:
                r[i+j] += pi * qj
    return tuple(r)

polymul_karatsuba(p, q)

Return the product of two polynomials.

\[ p q \]

Uses the Karatsuba algorithm.

Both arguments must be sequences.

TODO: complexity

See also
References
Source code in poly\standard\arithmetic.py 
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def polymul_karatsuba(p, q):
    """Return the product of two polynomials.

    $$
        p q
    $$

    Uses the Karatsuba algorithm.

    Both arguments must be sequences.

    TODO: complexity

    See also
    --------
    - for any implementation: [`polymul`][poly.standard.arithmetic.polymul]
    - other implementations:
    [`polymul_naive`][poly.standard.arithmetic.polymul_naive]

    References
    ----------
    - [Wikipedia - Karatsuba algorithm](https://en.wikipedia.org/wiki/Karatsuba_algorithm)
    """
    if not p or not q:
        return ()

    p, q = sorted((q, p), key=len) #p shorter, q longer
    ql, qu = q[:len(p)], q[len(p):]
    rl, ru = _polymul_karatsuba(p, ql), polymul_naive(p, qu)
    #return vecadd(rl, (0,)*len(p)+ru)
    return rl[:len(p)] + polyadd(rl[len(p):], ru)

polymulx(p, n=1, zero=0)

Return the product of polynomial p and a monomial of degree n.

\[ px^n \]

More efficient than polymul(p, polymonom(n)).

Complexity

There are no scalar arithmetic operations.

See also
References
Source code in poly\standard\arithmetic.py 
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
def polymulx(p, n=1, zero=0):
    """Return the product of polynomial `p` and a monomial of degree `n`.

    $$
        px^n
    $$

    More efficient than `polymul(p, polymonom(n))`.

    Complexity
    ----------
    There are no scalar arithmetic operations.

    See also
    --------
    - for polynomial factor: [`polymul`][poly.standard.arithmetic.polymul]
    - wraps: [`vector.vecrshift`](https://goessl.github.io/vector/functional/#vector.functional.utility.vecrshift)

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polymulx`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polymulx.html)
    """
    return vecrshift(p, n, zero=zero)

polypow(p, n, method='naive')

Return the polynomial p raised to the nonnegative n-th power.

\[ p^n \]

p must be a sequence.

Available methods are

TODO: mod parameter

See also
References
Source code in poly\standard\arithmetic.py 
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
def polypow(p, n, method='naive'):
    """Return the polynomial `p` raised to the nonnegative `n`-th power.

    $$
        p^n
    $$

    `p` must be a sequence.

    Available methods are 

    - [`naive`][poly.standard.polypow_naive] &
    - [`binary`][poly.standard.polypow_binary].

    TODO: mod parameter

    See also
    --------
    - implementations: [`polypow_naive`][poly.standard.arithmetic.polypow_naive],
    [`polypow_binary`][poly.standard.polypow_binary]
    - for sequence of powers: [`polypows`][poly.standard.arithmetic.polypows]

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polypow`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polypow.html)
    """
    match method:
        case 'naive':
            return polypow_naive(p, n)
        case 'binary':
            return polypow_binary(p, n)
        case _:
            raise ValueError('Invalid method')

polypow_naive(p, n, one=1)

Return the polynomial p raised to the nonnegative n-th power.

\[ p^n \]

Uses repeated multiplication.

p must be a sequence.

Complexity

For a polynomial of degree \(n\) and exponent \(k\) there will be

  • \(\begin{cases}\frac{n^2k(k-1)}{2}&n\ge0\\0&n\leq0\end{cases}\) scalar additions (add) &
  • \(\begin{cases}\frac{(nk+2)(n+1)(k-1)}{2}&k>0\\0&k=0\end{cases}\) scalar multiplications (mul).
See also
Source code in poly\standard\arithmetic.py 
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
def polypow_naive(p, n, one=1):
    r"""Return the polynomial `p` raised to the nonnegative `n`-th power.

    $$
        p^n
    $$

    Uses repeated multiplication.

    `p` must be a sequence.

    Complexity
    ----------
    For a polynomial of degree $n$ and exponent $k$ there will be

    - $\begin{cases}\frac{n^2k(k-1)}{2}&n\ge0\\0&n\leq0\end{cases}$ scalar additions (`add`) &
    - $\begin{cases}\frac{(nk+2)(n+1)(k-1)}{2}&k>0\\0&k=0\end{cases}$ scalar multiplications (`mul`).

    See also
    --------
    - for any implementation: [`polypow`][poly.standard.arithmetic.polypow]
    - other implementations:
    [`polypow_binary`][poly.standard.arithmetic.polypow_binary]
    - uses: [`polypows`][poly.standard.arithmetic.polypows]
    """
    return next(polypows(p, start=n, one=one))

polypow_binary(p, n, one=1)

Return the polynomial p raised to the nonnegative n-th power.

\[ p^n \]

Uses exponentiation by squaring.

p must be a sequence.

Complexity

For a polynomial of degree \(n\) and exponent \(k\) let

\[ \begin{aligned} k &= \sum_jb_j2^j \\ L &= \operatorname{bitlength}(k) \\ w &= \operatorname{popcount}(k) \end{aligned} \]

where \(j_0\) is the least significant \(1\)-bit position of \(k\), \(\operatorname{bitlength}(k)\) is the number of bits of the binary representation of \(k\) and \(\operatorname{popcount}(k)\) is the number of \(1\)-bits.

Further define

\[ \begin{aligned} A(k) &= \frac{4^L-1}{3}+\sum_{\substack{i<j \\ b_i=b_j=1}}2^{i+j} \\ B(k) &= 2(2^L-1)+k-2^{j_0}+\sum_{j\mid b_j=1}2^j\operatorname{popcount}(k\gg(j+1)) \\ C(k) &= L+w-1. \end{aligned} \]
L = k.bit_length()
w = k.bit_count()
C = L + w - 1
B = 2*(2**L-1) + k-2**((k&-k).bit_length()-1) + sum(2**j * (k>>(j+1)).bit_count() for j in range(L) if ((k>>j)&0x01)==0x01)
A = (4**L-1)//3 + sum(2**(i+j) for j in range(L) for i in range(j) if ((k>>j)&0x01)==((k>>i)&0x01)==0x01)

Then there will be

  • \(A(k)n^2\) scalar additions (add) &
  • \(A(k)n^2+B(k)n+C(k)\) scalar multiplications (mul).
See also
References
Source code in poly\standard\arithmetic.py 
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
def polypow_binary(p, n, one=1):
    r"""Return the polynomial `p` raised to the nonnegative `n`-th power.

    $$
        p^n
    $$

    Uses exponentiation by squaring.

    `p` must be a sequence.

    Complexity
    ----------
    For a polynomial of degree $n$ and exponent $k$ let

    $$
        \begin{aligned}
            k &= \sum_jb_j2^j \\
            L &= \operatorname{bitlength}(k) \\
            w &= \operatorname{popcount}(k)
        \end{aligned}
    $$

    where $j_0$ is the least significant $1$-bit position of $k$,
    $\operatorname{bitlength}(k)$ is the number of bits of the binary
    representation of $k$ and $\operatorname{popcount}(k)$ is the number of
    $1$-bits.

    Further define

    $$
        \begin{aligned}
            A(k) &= \frac{4^L-1}{3}+\sum_{\substack{i<j \\ b_i=b_j=1}}2^{i+j} \\
            B(k) &= 2(2^L-1)+k-2^{j_0}+\sum_{j\mid b_j=1}2^j\operatorname{popcount}(k\gg(j+1)) \\
            C(k) &= L+w-1.
        \end{aligned}
    $$

    ```python
    L = k.bit_length()
    w = k.bit_count()
    C = L + w - 1
    B = 2*(2**L-1) + k-2**((k&-k).bit_length()-1) + sum(2**j * (k>>(j+1)).bit_count() for j in range(L) if ((k>>j)&0x01)==0x01)
    A = (4**L-1)//3 + sum(2**(i+j) for j in range(L) for i in range(j) if ((k>>j)&0x01)==((k>>i)&0x01)==0x01)
    ```

    Then there will be

    - $A(k)n^2$ scalar additions (`add`) &
    - $A(k)n^2+B(k)n+C(k)$ scalar multiplications (`mul`).

    See also
    --------
    - for any implementation: [`polypow`][poly.standard.arithmetic.polypow]
    - other implementations:
    [`polypow_naive`][poly.standard.arithmetic.polypow_naive]

    References
    ----------
    - [Wikipedia - Exponentiation by squaring](https://en.wikipedia.org/wiki/Exponentiation_by_squaring)
    - Sequence $C(k)$: [A056791](https://oeis.org/A056791)
    """
    if n == 0:
        return (one,)
    r = None
    while n:
        if n % 2 == 1:
            r = polymul_naive(r, p) if r is not None else p
        p = polymul_naive(p, p)
        n //= 2
    return r

polypows(p, start=0, one=1)

Yield the powers of the polynomial p.

\[ (p^n)_{n\in\mathbb{N}_0} = (1, p, p^2, \dots) \]

Uses iterative multiplication to calculate powers consecutively.

p must be a sequence.

Notes

Was first .evaluation.polycoms in analogy to polyvals but then the submodules arithmetic & evaluation would not be separable (.evaluation.polycoms uses .arithmetic.polymul and .arithmetic.polypow uses .evaluation.polycoms).

See also
Source code in poly\standard\arithmetic.py 
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
def polypows(p, start=0, one=1):
    r"""Yield the powers of the polynomial `p`.

    $$
        (p^n)_{n\in\mathbb{N}_0} = (1, p, p^2, \dots)
    $$

    Uses iterative multiplication to calculate powers consecutively.

    `p` must be a sequence.

    Notes
    -----
    Was first `.evaluation.polycoms` in analogy to
    [`polyvals`][poly.standard.polyvals] but then the submodules `arithmetic` &
    `evaluation` would not be separable (`.evaluation.polycoms` uses
    `.arithmetic.polymul` and `.arithmetic.polypow` uses
    `.evaluation.polycoms`).

    See also
    --------
    - used by: [`polypow_naive`][poly.standard.arithmetic.polypow_naive], [`polycom_iterative`][poly.standard.polycom_iterative]
    - for scalar arguments: [`polyvals`][poly.standard.polyvals]
    """
    if start <= 0:
        yield (one,)
    yield (q := reduce_default(polymul_naive, repeat(p, start), initial=MISSING, default=p))
    while True:
        q = polymul_naive(q, p)
        yield q

calculus

polyder(p, k=1)

Return the k-th derivative of polynomial p.

\[ p^{(k)} \]
Complexity

For the \(k\)-th derivative of a polynomial of degree \(n\) there will be

  • \(\begin{cases}n-k+1&k\le n\\0&k>n\end{cases}\) scalar multiplications with integers (rmul).
Notes

For monomials:

\[ \begin{aligned} \frac{d}{dx}x^n &= nx^{n-1} \\ \frac{d^2}{dx^2}x^n &= n(n-1)x^{n-2} \\ &\vdots \\ \frac{d^k}{dx^k}x^n &= n(n-1)\cdots(n-k+1)x^{n-k} \\ &= (n)_kx^{n-k} \\ &= \frac{n!}{(n-k)!}x^{n-k} \\ &= {}_nP_kx^{n-k} \end{aligned} \]

And for polynomials:

\[ \begin{aligned} \frac{d^k}{dx^k}p(x) &= \frac{d^k}{dx^k}\sum_{l=0}^np_lx^l \\ &= \sum_{l=k}^np_l\,{}_lP_kx^{l-k} \\ &= \sum_{l=0}^{n-k}p_{l+k}\,{}_{l+k}P_kx^l \end{aligned} \]

Where \((n)_k\) is the Falling factorial and \({}_nP_k\) is the number of k-permutations of n with \((n)_k=\frac{n!}{(n-k)!}={}_nP_k\) (falling factorials are used because their definition appears in the derivation; permutations are used because a fast implementation is provided by math.perm).

References
Source code in poly\standard\calculus.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def polyder(p, k=1):
    r"""Return the `k`-th derivative of polynomial `p`.

    $$
        p^{(k)}
    $$

    Complexity
    ----------
    For the $k$-th derivative of a polynomial of degree $n$ there will be

    - $\begin{cases}n-k+1&k\le n\\0&k>n\end{cases}$ scalar multiplications with integers (`rmul`).

    Notes
    -----
    For monomials:

    $$
        \begin{aligned}
            \frac{d}{dx}x^n &= nx^{n-1} \\
            \frac{d^2}{dx^2}x^n &= n(n-1)x^{n-2} \\
            &\vdots \\
            \frac{d^k}{dx^k}x^n &= n(n-1)\cdots(n-k+1)x^{n-k} \\
            &= (n)_kx^{n-k} \\
            &= \frac{n!}{(n-k)!}x^{n-k} \\
            &= {}_nP_kx^{n-k}
        \end{aligned}
    $$

    And for polynomials:

    $$
        \begin{aligned}
            \frac{d^k}{dx^k}p(x) &= \frac{d^k}{dx^k}\sum_{l=0}^np_lx^l \\
            &= \sum_{l=k}^np_l\,{}_lP_kx^{l-k} \\
            &= \sum_{l=0}^{n-k}p_{l+k}\,{}_{l+k}P_kx^l
        \end{aligned}
    $$

    Where $(n)_k$ is the [Falling factorial](https://en.wikipedia.org/wiki/Falling_and_rising_factorials#Properties)
    and ${}_nP_k$ is the [number of k-permutations of n](https://en.wikipedia.org/wiki/Permutation#k-permutations_of_n)
    with $(n)_k=\frac{n!}{(n-k)!}={}_nP_k$ (falling factorials are
    used because their definition appears in the derivation; permutations are
    used because a fast implementation is provided by `math.perm`).

    References
    ----------
    - [more_itertools.polynomial_derivative](https://more-itertools.readthedocs.io/en/stable/api.html#more_itertools.polynomial_derivative)
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polyder`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyder.html)
    - [Wikipedia - Potenzregel - Höhere Ableitung einer Potenzfunktion mit natürlichem Exponenten](https://de.wikipedia.org/wiki/Potenzregel#H%C3%B6here_Ableitung_einer_Potenzfunktion_mit_nat%C3%BCrlichem_Exponenten)
    """
    return vechadamard((perm(l+k, k) for l in count()), islice(p, k, None))

polyantider(p, c=0, b=0)

Return the antiderivative of polynomial p.

\[ \int_b^xp(y)\,\mathrm{d}y+c \]

TODO: Higher antiderivatives, complexity

Notes

Let

\[ I_{b, c}[f](x) = \int_b^xf(y)\,\mathrm{d}y+c. \]

Then we have for monomials:

\[ \begin{aligned} I_{b,c}[\cdot^n](x) &= \int_b^xy^n\,\mathrm{d}y+c \\ &= \left.\frac{y^{n+1}}{n+1}\right|_{y=b}^x+c \\ &= \frac{x^{n+1}}{n+1}-\frac{b^{n+1}}{n+1}+c \end{aligned} \]

For polynomials:

\[ \begin{aligned} I_{b,c}[p](x) &= I_{b,c}\left[\sum_ka_k\cdot^k\right](x) \\ &= \int_b^x\sum_ka_ky^k\,\mathrm{d}y+c \\ &= \sum_ka_k\left(\frac{x^{k+1}}{k+1}-\frac{b^{k+1}}{k+1}\right)+c \\ &= \sum_ka_{k-1}\left(\frac{x^k}{k}-\frac{b^k}{k}\right)+c \\ &= \sum_k\frac{a_{k-1}}{k}x^k-\sum_k\frac{a_{k-1}}{k}b^k+c \end{aligned} \]
Notes

Integration is called antiderivative (antider) instead of integrate (int) to avoid keyword collisions.

References
Source code in poly\standard\calculus.py 
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def polyantider(p, c=0, b=0):
    r"""Return the antiderivative of polynomial `p`.

    $$
        \int_b^xp(y)\,\mathrm{d}y+c
    $$

    TODO: Higher antiderivatives, complexity

    Notes
    -----
    Let

    $$
        I_{b, c}[f](x) = \int_b^xf(y)\,\mathrm{d}y+c.
    $$

    Then we have for monomials:

    $$
        \begin{aligned}
            I_{b,c}[\cdot^n](x) &= \int_b^xy^n\,\mathrm{d}y+c \\
            &= \left.\frac{y^{n+1}}{n+1}\right|_{y=b}^x+c \\
            &= \frac{x^{n+1}}{n+1}-\frac{b^{n+1}}{n+1}+c
        \end{aligned}
    $$

    For polynomials:

    $$
        \begin{aligned}
            I_{b,c}[p](x) &= I_{b,c}\left[\sum_ka_k\cdot^k\right](x) \\
            &= \int_b^x\sum_ka_ky^k\,\mathrm{d}y+c \\
            &= \sum_ka_k\left(\frac{x^{k+1}}{k+1}-\frac{b^{k+1}}{k+1}\right)+c \\
            &= \sum_ka_{k-1}\left(\frac{x^k}{k}-\frac{b^k}{k}\right)+c \\
            &= \sum_k\frac{a_{k-1}}{k}x^k-\sum_k\frac{a_{k-1}}{k}b^k+c
        \end{aligned}
    $$

    Notes
    -----
    Integration is called antiderivative (`antider`)
    instead of integrate (`int`) to avoid keyword collisions.

    References
    ----------
    - `numpy` equivalent: [`numpy.polynomial.polynomial.polyint`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyint.html)
    """
    P = vechadamardtruediv(p, count(1))
    return vecrshift(P, 1, zero=c-b*polyval_iterative(P, b) if b else c)

conversion

polysympify(p, x=x)

Return the coefficient iterable p as a sympy.Poly.

Source code in poly\standard\conversion.py 
10
11
12
def polysympify(p, x=x):
    """Return the coefficient iterable `p` as a [`sympy.Poly`](https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly)."""
    return sp.Poly.from_list(tuple(reversed(tuple(p))), x)

polyunsympify(p)

Return sympy.Poly(p) as a coefficient tuple.

Source code in poly\standard\conversion.py 
14
15
16
def polyunsympify(p):
    """Return [`sympy.Poly(p)`](https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly) as a coefficient `tuple`."""
    return tuple(reversed(p.all_coeffs()))