Symbolic implementation of a polynomial using numpy and bare python.

nppoly.py

import numpy as np
from functools import reduce, lru_cache


@lru_cache(maxsize=None)
def binom(n, k):
    if k == 0: return 1
    if n == k: return 1
    return binom(n - 1, k - 1) + binom(n - 1, k)


class Polynomial:
    """
    Implementes sum_j theta_j x^j mathematical object.
    """
    def __init__(self, *coeffs):
        c = np.trim_zeros(np.array(coeffs), trim='f')
        self._coeffs = c if c.shape[0] else np.zeros(1)


    @classmethod
    def monomial(cls, power, value):
        """
        Defines value * x^power mathematical object.
        """
        coeffs = [0] * (power + 1)
        coeffs[0] = value
        return cls(*coeffs)


    def __str__(self):
        c = np.linspace(len(self), 0, num=len(self) + 1, dtype='int')
        return np.array2string(self._coeffs, separator='x{} ') \
                 .format(*c).strip('[]')


    def __repr__(self):
        return f'Polynomial{tuple(self._coeffs)}'


    @property
    def shape(self):
        return self._coeffs.shape[0]
    

    def __len__(self):
        return self.shape - 1


    def __call__(self, x):
       x = x.reshape(-1, 1)
       exponents = np.linspace(len(self), 0, num=self.shape)
       return (x ** exponents * self._coeffs).sum(axis=1)


    def __add__(self, other):
        if isinstance(other, (int, float)):
            result = self._coeffs + np.append(np.zeros(len(self)), other)
            return Polynomial(*tuple(result))
        elif isinstance(other, Polynomial):
            result = np.zeros(max(self.shape, other.shape))
            result[len(result) - self.shape:] = self._coeffs
            result[len(result) - other.shape:] += other._coeffs
            return Polynomial(*tuple(result))
        else:
            raise(NotImplementedError)


    def __radd__(self, other):
        return self.__add__(other)


    def __sub__(self, other):
        if isinstance(other, (int, float)):
            result = self._coeffs - np.append(np.zeros(len(self)), other)
            return Polynomial(*tuple(result))
        elif isinstance(other, Polynomial):
            result = np.zeros(max(self.shape, other.shape))
            result[len(result) - self.shape:] = self._coeffs
            result[len(result) - other.shape:] -= other._coeffs
            return Polynomial(*tuple(result))
        else:
            raise(NotImplementedError)


    def __mul__(self, other):
        if isinstance(other, (int, float)):
            coeffs = self._coeffs * other
            return Polynomial(*coeffs)
        elif isinstance(other, Polynomial):
            C = np.zeros((other.shape, self.shape + other.shape - 1))
            for i in range(other.shape):
                C[i].put(range(i, len(self._coeffs) + i), 
                         self._coeffs * other._coeffs[i])
            return Polynomial(*tuple(C.sum(axis=0)))
        else:
            raise(NotImplementedError)
            
            
    def __rmul__(self, other):
        return self.__mul__(other)


    def __matmul__(self, other):
        return self.__mul__(other)


    def __rsub__(self, other):
        return self.__sub__(other)


    def __pow__(self, n):
        if isinstance(n, int):
            if n == 0:
                return Polynomial(1)
            elif n > 0:
                return reduce(lambda x, y: x * y, [self] * n) 
            else:
                raise(ValueError)
        else:
            raise(NotImplementedError)

    
    def __truediv__(self, other):
        coeffs = self._coeffs / other
        return Polynomial(*tuple(coeffs))


    def __floordiv__(self, other):
        coeffs = self._coeffs // other
        return Polynomial(*tuple(coeffs))


    def __lshift__(self, value):
        coeffs = np.append(self._coeffs, np.zeros(value))
        return Polynomial(*tuple(coeffs))

    
    def __rshift__(self, value):
        coeffs = self._coeffs[slice(None, -value)]
        return Polynomial(*tuple(coeffs))

    
    def __eq__(self, other):
        if not isinstance(other, Polynomial):
            raise(NotImplementedError)
        return tuple(self._coeffs) == tuple(other._coeffs)


    def __getitem__(self, *powers):
        p = len(self) - np.array(powers)
        c = np.zeros(self.shape)
        c.put(p, np.ones_like(p))
        return Polynomial(*tuple(c * self._coeffs))


    def __setitem__(self, powers, values):
        p = len(self) - np.array(powers)
        self._coeffs.put(p, values)

    
    def __delitem__(self, *powers):
        p = len(self) - np.array(powers)
        self._coeffs.put(p, np.zeros_like(p))
        self.__init__(*tuple(self._coeffs))


    def grad(self, value=None):
        coeffs = (self >> 1)._coeffs * np.linspace(len(self), 1, num=len(self))
        if value is None:
            return Polynomial(*tuple(coeffs))
        return Polynomial(*tuple(coeffs))(value)