Polynomial methods and operators

Coercion into polynomial rings

First let’s deal with the rather interesting task of overloading opCall for our poly class. Here is the method:

static poly!(T, X) opCall(S)(S p)
{
   static if (__traits(hasMember, S, "coeff"))
   {
      static if(is(typeof(p.coeff) == T[]))
         return p;
      else
      {
         auto c = T(p);
         return new poly!(T, X)(1);
      }
   } else
      return new poly!(T, X)([p]);
}

There are two broad cases we are dealing with here. The first is if we are coercing a polynomial over a particular ring to a polynomial over T. The other case is if we are coercing a value (such as an int or ZZ) into the polynomial ring over T. The latter case is easy to deal with, as we already have a constructor which takes an array of coefficients and turns it into a polynomial, so we simply call that.

To detect if we are in the first case, we need to determine if p is a polynomial. If it is, it will have a “coeff” field. This can be detected using the hasMember trait. Again, this is a static compile time decision and so can go in a “static if”.

If we do have a polynomial, there are again two cases. Either that polynomial is already a polynomial over T, in which case no coercion is required and we can simply return p.

In the alternative case, we must be coercing from an element of a smaller ring. Thus, we should consider p to be a constant coefficient in the new ring over T, i.e. we first need to coerce p into the ring T, then create a polynomial whose only coefficient is this value.

The above code can probably be made more succinct by combining cases, but what we have is comprehensible enough.

Overloading addition for polynomials

The benchmark that we give at the end of this blog series will only require addition of polynomials. For that we overload opAdd. It’s a very straightforward polynomial addition function.

There are two cases: either we are adding two polynomials over T or we are adding something else to a polynomial over T. In the latter case we first coerce to a polynomial over T using the coercion we just implemented. The former case is a very straightforward polynomial addition function.

poly!(T, X) opAdd(S)(S b)
{
   static if (is(S == poly!(T, X)))
   {
      size_t i;
      auto r = new poly!(T, X)();

      r.coeff.length = max(coeff.length, b.coeff.length);

      for (i = 0; i < min(coeff.length, b.coeff.length); i++)
         r.coeff[i] = coeff[i] + b.coeff[i];

      for ( ; i < coeff.length; i++)
         r.coeff[i] = coeff[i];

      for ( ; i < b.coeff.length; i++)
         r.coeff[i] = b.coeff[i];

      r.normalise();

      return r;
   } else
   {
      auto t = poly!(T, X)(b);
      return this + t;
   }
}

The D standard library provides the functions min and max in the module std.algorithm.

Addition proceeds exactly how one might expect it to. A new polynomial is created to hold the output, then the input coefficients are added, being careful in the case that one of the polynomials is longer than the other.

Normalisation

It may also be that addition of coefficients causes cancellation. We prefer the top coefficient of our polynomial to be nonzero, so we normalise it so that this is the case (the zero polynomial simply has length zero).

Here is the normalise method:

void normalise()
{
   while (coeff.length && coeff[$ - 1] == 0)
      coeff.length = coeff.length - 1;
}

D allows the coefficient array to be reallocated simply by adjusting its length with an assignment, which is a neat trick. Also the length of an array is given by the $ operator inside the square brackets, which saves typing.

Equality operator

Another thing which our code assumes so far is that we have a method to compare a polynomial with zero. For this, we need to overload the relevant comparison operator opEquals in the case of an int parameter. The only tricky part here is that there are two cases: when we wish to compare with zero (in which case we want to test for a polynomial of length zero) and the case where we wish to compare with a nonzero value.

bool opEquals(int p)
{
   if (p == 0) return coeff.length == 0;
   else return coeff.length == 1 && coeff[0] == p;
}

The toString method

Finally, we need to be able to print polynomials, which means we need to override the toString method.

override string toString()
{
   string s = "";
   size_t i;

   if (coeff.length)
      for (i = coeff.length - 1; i > 1; i--)
         s = s ~ br(i, "(") ~ coeff[i].toString() ~ br(i, ")")
              ~ "*" ~ X.stringof[1..$-1] ~ "^" ~ to!(string)(i) ~ " + ";

   if (coeff.length > 1)
      s = s ~ br(1, "(") ~ coeff[1].toString() ~ br(1, ")")
            ~ "*" ~ X.stringof[1..$-1] ~ " + ";

   if (coeff.length)
      s = s ~ br(0, "(") ~ coeff[0].toString() ~ br(0, ")");
   else
      s = "0";

   return s;
}

Here we make use of a really neat trick. We need to be able to print the variable for the polynomial. But recall that we placed this information in the type, in the parameter X. To print this we make use of D’s stringof method.

Actually, stringof is a property method, which means that it can be accessed like a field of X, even though it is a method. All expressions and types in D have such a stringof property method (users can write their own as well). This will give us the string “x”, including the double quotation marks, which is a bit odd, so we shave these off by taking a slice of the string.

Another thing we need to do is to put parentheses around any coefficients which are themselves polynomials, otherwise (x^2 + 1)*y would be printed as x^2 + 1*y, for example.

To handle this we implement a function br which checks if parentheses are needed for a given coefficient:

string br(size_t i, string s)
{
   static if (is(typeof(coeff[i].coeff)))
      if (coeff[i].coeff.length > 1)
         return s;

   return "";
}

We don’t print parentheses unless the coefficient is a polynomial and of length greater than one.

Note on performance

At this point we actually have everything we need to implement a multivariate polynomial addition benchmark. It’s very important to realise this. All the code that comes after this point is included for convenience of the user and does not increase the runtime of our benchmark in any significant way.

It is by no means obvious upon first viewing of the code that this is the case. Indeed it looks like the code that appears in the following posts is used every time a polynomial is created. This is true, but only when the user creates a polynomial, and that isn’t going to happen often because fingers!!

This all becomes clearer if we realise that we can already create polynomials using the constructors we have, and add those polynomials together in a loop, etc.

The code that is described from this point on is only used to build types and polynomials, once and for all, in a more convenient manner than with the existing constructors!

‹‹‹ Previous Post Next Post ›››

Leave a Reply