Page MenuHomec4science
Files F62985789

polynomial_algebra.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 23:14

polynomial_algebra.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import Np2D, Tuple
from .vander import polyvander
from rrompy.utilities.numerical import pseudoInverse
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['changePolyBasis', 'polyTimes', 'polyDivide']
def changePolyBasis(P:Np2D, dim : int = None, basis0 : str = "MONOMIAL",
basisF : str = "MONOMIAL") -> Np2D:
if basis0 == basisF: return P
if dim is None: dim = P.ndim
if basis0 != "MONOMIAL" and basisF != "MONOMIAL":
return changePolyBasis(changePolyBasis(P, dim, basis0, "MONOMIAL"),
dim, "MONOMIAL", basisF)
basisD = basisF if basis0 == "MONOMIAL" else basis0
R = copy(P)
N = np.max(P.shape[: dim]) - 1
vander = polyvander([0], N, basisD, [list(range(N + 1))])
if basis0 == "MONOMIAL": vander = pseudoInverse(vander)
for j in range(dim):
R = np.tensordot(vander, R, (-1, j))
return R
def polyTimes(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
Qbasis : str = "MONOMIAL", Rbasis : str = "MONOMIAL") -> Np2D:
if not isinstance(P, (np.ndarray,)): P = np.array(P)
if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
if dim is None: dim = P.ndim
if dim <= 0: return
R = np.zeros([x + y - 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
dtype = P.dtype)
if dim == 1:
for j, Qj in enumerate(Q):
R[j : j + len(P)] = R[j : j + len(P)] + Qj * P
else:
for j, Qj in enumerate(Q):
for l, Pl in enumerate(P):
R[j + l] = R[j + l] + polyTimes(Pl, Qj, dim - 1)
return changePolyBasis(R, dim, "MONOMIAL", Rbasis)
def polyDivide(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
Qbasis : str = "MONOMIAL",
Rbasis : str = "MONOMIAL") -> Tuple[Np2D, Np2D]:
if not isinstance(P, (np.ndarray,)): P = np.array(P)
if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
Pc = copy(P)
Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
if dim is None: dim = P.ndim
if dim <= 0: return
R = np.zeros([x - y + 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
dtype = P.dtype)
if dim == 1:
for i in range(len(R) - 1, -1, -1):
R[i] = Pc[-1] / Q[-1]
Pc = Pc[: -1]
for j, Qj in enumerate(Q[::-1]):
if j > 0: Pc[-j] = Pc[-j] - R[i] * Qj
else:
raise RROMPyException(("Quotient of multivariate polynomials not "
"supported."))
return (changePolyBasis(R, dim, "MONOMIAL", Rbasis),
changePolyBasis(Pc, dim, "MONOMIAL", Rbasis))

Event Timeline

c4science · Help