An internal DSL for polynomial rings

Syntactic sugar

At present there are three problems for the user if they wish to create polynomials.

  • Annotating the types of the polynomials can be a real chore, e.g. poly!(poly!(poly!(ZZ, “x”), “y”), “z”)
  • We can’t yet build polynomials up using algebraic notation f = (3*x^3 + 2*x + 1)*y + (2*x + 1), etc.

What we would really like is to be able to give a name, R say, to poly!(ZZ, “x”) and then build S say, as poly!(R, “y”). This would also make coercion easier to annotate, as we could write S(f) instead of poly!(poly!(ZZ, “x”), “y”)(f).

Type aliases

Fortunately, D has a means of assigning an alias to an existing type. Here is how we can use it:

alias poly!(ZZ, "x") R;
alias poly!(R, "y") S;

This solves the coercion and type annotation problem.

Polynomial rings and generators

The second thing we would like to do boils down to being able to define x to be the variable in our ring R, i.e. the generator of our ring. The logical way to do this is to have a function/macro, PolynomialRing say, which takes a ring T and a variable string X and outputs a tuple containing a type R representing the ring and a value x representing the generator of this ring. Then we could simply type:

auto (R, x) = PolynomialRing(ZZ, "x");

Quite surprisingly, D does not allow tuples to be returned from functions and does not seem to have destructuring binds. Therefore, I don’t see a way to do this.

Even if D offered this feature, there would be no requirement that it provide the ability to return a type from a function and then use that type in later code. This is something you expect in dynamic languages. Taking type information over to the run time side and back to the compile side is not something you expect to be supported in a statically typed language.

However, what can be offered instead is the ability to define macros. These allow the equivalent of a function to be defined which is expanded out at compile time. In combination with tuples and destructuring binds, many new things become possible. Sadly, I couldn’t see how to do this in D.

Anyhow, instead or returning a tuple, I chose to implement a method gen() of the poly class which returns something representing the generator x of the polynomial ring. One way to do this is to have gen() simply return the polynomial poly!(T, X)([0, 1]).

A monomial class

However, this becomes quite complex without implementing general multiplicaton of polynomials. So to make life easier for myself I decided to add an additional class called monomial which allows a power of the generator of a polynomial ring to be represented. For example, x^7 would be represented as a monomial!(T, X) containing the value 7.

This may seem like a roundabout way of doing things, but it makes implementation easier.

Doing things this way can also have a performance benefit when constructing very long polynomials. For example, suppose we wanted to construct the polynomial 3*x^1000000 + 2*x^100000. If the monomial type contained a power of x and a coefficient, we could define addition of two monomials in an efficient way which does not require creation of two extremely long arrays and another long array representing their sum.

I decided not to include the coefficient in the monomial type, so my monomials only contain the exponent. This just made the implementation easier and doesn’t affect the benchmark in any way as we only use short polynomials.

Here is the monomial class:

class monomial(T, string X)
   size_t exp;

   this(size_t p) { exp = p; }

The constructor takes a size_t representing the exponent, e.g. x^7 would require exp to be set to 7.

Exponentiation operator

The first thing we want to be able to do with a monomial is raise it to a power. Initially x will be the monomial with exp = 1, thus to create x^7 we simply want to multiply that value of exp by 7.

monomial!(T, X) opPow(size_t p)
   return new monomial!(T, X)(p*exp);

At first I made the mistake of overloading the ^ operator, called opXor. However, I soon remembered that this has left-to-right associativity and a lower prededence than opMul. I needed opPow, which is denoted ^^ in D. This has the right associativity and precedence for exponentiation. Some languages such as Python use ** instead of ^^ for exponentiation. Sage uses ^ in interactive sessions because most mathematicians are so used to it. I understand that it makes use of a very thin preparser to allow this.

Multiplying a monomial and a coefficient

Now we need to think about how to create 3*x^7, or indeed (2*x + 1)*y^8, which in each case requires multiplication of a monomial by some coefficient. To accomplish this we coerce the coefficient into the ring T and turn the monomial into a polynomial with that coefficient in the right place in the coeff array.

poly!(T, X) opMul(S)(S c)
   auto t = new poly!(T, X)(this);
   t.coeff[exp] = T(c);
   return t;

Printing monomials

Finally we need to be able to print a monomial on its own.

override string toString()
   if (exp > 1)
      return X.stringof[1..$-1] ~ "^" ~ to!(string)(exp);
   else if (exp)
      return X.stringof[1..$-1];
      return "1";

The generator of the polynomial ring

Now we can add a gen() method to our poly class. It will create the monomial with exp = 1.

static monomial!(T, X) gen()
   return new monomial!(T, X)(1);

We make it a static method of the class so that we can write, for example:

alias poly!(ZZ, "x") R;

auto x = R.gen();

All that remains is to implement a poly constructor which accepts a monomial and overload the poly class’s opCall operator to accept a monomial.

These are simply extra cases in the “static if” statements for the relevant methods.

The code so far

All of the code so far can be found here:

We are finally ready to implement our very simple benchmark.

‹‹‹ Previous Post Next Post ›››

Leave a Reply