Home > OS >  Usage example of polynomial-algebra
Usage example of polynomial-algebra

Time:07-24

I'm looking for a Haskell library to manipulate multivariate polynomials. I found polynomial-algebra but there's no example and the doc is too complicated for me. I'd like some examples of:

  • constructing a multivariate polynomial
  • adding and multiplying two multivariate polynomials

with this library.

In addition, I'd like to know whether this library allows to deal with symbolic coefficients, e.g.

(1 a)*x1**2   x1*x2   a

where a is a literal variable (I'm able to do that with SymPy in Python and with AbstractAlgebra in Julia). That is to say the polynomial has coefficients in the field (or ring?) Q(a) or Z(a).

CodePudding user response:

Here are some things I can tell by reading the polynomial-algebra docs you linked:

  • You are offered a choice of several ways to describe the variables
  • The coefficient is abstracted over, and can be any type you like
  • The polynomial types all support Num so long as the coefficient type you choose is a Ring, so you add and multiply them using the ordinary ( ) and (*) functions.
  • There is built in support for using Int, Integer, and Rational as coefficients. The polynomial type itself also supports Ring, so you can use a polynomial in some different variable as a coefficient.

Given all that, it seems like this library meets your needs. Let's read the docs further and see how to actually use it.

But before we get to your example multivariate polynomial, let's start with a simple univariate example (which you'll need anyway to represent a later). First we need DataKinds for the type-level strings that denote the variables:

{-# LANGUAGE DataKinds #-}
import qualified Math.Algebra.Polynomial.Univariate as Uni
import Math.Algebra.Polynomial.Univariate (U(..))

import qualified Math.Algebra.Polynomial.FreeModule as FM
import Math.Algebra.Polynomial.Pretty (pretty)

Now, let's define y = 1 x^2 as our simple polynomial:

y :: Uni.Univariate Int "x"
y = Uni.Uni (FM.fromList [(U 0, 1), (U 2, 1)])

And write a main function that prints y and its square in a human-readable way:

main = do
  putStrLn $ (pretty y)    " squared equals"
  putStrLn (pretty (y * y))

Observe that we get the result we wanted:

1   x^2 squared equals
1   2*x^2   x^4

Now, on to multivariate polynomials. There are various ways you could choose to represent (1 a)*x1**2 x1*x2 a. The main decision is how to represent the various xs. What seems simplest to me is to use Math.Algebra.Polynomial.Multivariate.Infinite to support xn for an arbitrarily large n. For slightly more efficient operations (storing the monomials in an unboxed array instead of a list), you could instead use Math.Algebra.Polynomial.Multivariate.Indexed.

This time I'll just include the whole source file as one block:

{-# LANGUAGE DataKinds #-}

import qualified Math.Algebra.Polynomial.Univariate as Uni
import Math.Algebra.Polynomial.Univariate (U(..))

import qualified Math.Algebra.Polynomial.Multivariate.Infinite as Inf
import Math.Algebra.Polynomial.Multivariate.Infinite (XInf(..))

import qualified Math.Algebra.Polynomial.FreeModule as FM
import Math.Algebra.Polynomial.Pretty (pretty)

y :: Inf.Poly (Uni.Univariate Int "a") "x"
y = Inf.Poly
    (FM.fromList
     [ (XInf [], Uni.Uni (FM.singleton (U 1) 1)) -- a
     , (XInf [1, 1], Uni.Uni (FM.singleton (U 0) 1)) -- x1 * x2
     , (XInf [2], Uni.Uni (FM.fromList [(U 0, 1), (U 1, 1)])) -- (1   a)x1^2
     ])

main = do
  putStrLn $ (pretty y)    " squared equals"
  putStrLn (pretty (y * y))

We get a result that looks plausible to me, but I haven't done the expansion by hand to check:

a*(1)   1*x1*x2   1   a*x1^2 squared equals
a^2*(1)   2*a*x1*x2   2*a   2*a^2*x1^2   1*x1^2*x2^2   2   2*a*x1^3*x2   1   2*a   a^2*x1^4

One issue you can see with this approach is that the pretty-printing doesn't include parens around the coefficients, so it prints 1 a*x1^2 instead of (1 a)x1^2. I'm not sure of the best way to address this. It seems like there ought to be some place where you can insert a newtype around the coefficient type that causes it to always render with parens around it, but I couldn't find it.

  • Related