The Polynomial Object of the Computer Algebra Kit

David Stes

March 6, 1996

Abstract:

Brad Cox' book on OBJECTIVE-C [Cox, 1986] is a plea for building software components that look like Smalltalk classes, but that are implemented in C instead of Smalltalk. This concept has proven to be useful for writing UNIX filters, such as compilers or text formatting tools, or for developing software with sophisticated, graphical user interfaces. A new application of the technology can be found in another field of notoriously big and complex software systems : the construction of interactive software with symbolical user interfaces. The kernels of these systems can be composed out of Smalltalk-like objects, while the components themselves can be written in C. We[*] have implemented, over the last three years[*], a set of -- essentially -- as many OBJECTIVE-C objects: Integer , Polynomial and Matrix , and we would like to comment in this paper on the design of the single most important object, the Polynomial object of the Computer Algebra Kit library.

Smalltalk

The phrase ``components that look like Smalltalk classes'' deserves some further elaboration. Smalltalk programming style [Goldberg and Robson, 1983] is characterized by the use of such notions as objects, messages, blocks, collections and streams. It suffices to say that, in Smalltalk, everything's an Object, and Collections or Sequences are groups of arbitrary objects. For example, Sequences are groups of objects meant to be sequentially accessed.

Using the OBJECTIVE-C syntax, some hybrid of C and Smalltalk :

aSequence = [aPolynomial eachTerm];
This statement creates -- given a Polynomial object -- a new Sequence object that consists of the Term objects of the polynomial; a term consists of a coefficient multiplied by a symbol raised to some exponent.

The following code snippet iterates over the terms in the sequence, and computes the sum of the squares of the coefficients of the terms :

while (aTerm = [aSequence nextMember]) {
	sqr = [[aTerm coefficient] square];
	sum = [sum add:sqr];
}

The coefficients and symbols of the terms, that make up a polynomial, are arbitrary objects, said to be loosely coupled to the Polynomial object : it's decided only at run-time what kind of coefficients (integer, fractional, modular, floating-point etc.) the polynomial actually has (or what kind of symbols), and it is only at run-time that it is decided what the squaring or addition of these coefficients might look like.

Expression Trees

Software systems that interact with the user through the use of symbolic expressions, generally represent user input in a tree that is in a one-to-one relationship with the input. Different systems use different notations, and the users of these systems find some particular notation most convenient for their application. Typically, an expression tree consists of integer or floating-point nodes, of symbolic parameters, of the action commands specific to the system and so on.

Computer algebra and polynomial arithmetic is just one component in the construction of such software packages with symbolic user interfaces : in many applications, the use of symbolic expressions is the glue between the more essential, computational, numerical or graphical components of the system. How polynomials enter the picture, can be seen when considering, for example, the processing of the statement :

Simplification towards :

can be understood as computing with a polynomial, with integer coefficients, in the symbols and .

For the Polynomial object to be reusable between the kernels of different systems -- using different expression trees -- or simply for sharing the code between different versions of a same system, it's important that the Polynomial object is not tightly coupled to a particular system : in addition to the fact that the object must be capable of working with different kinds of coefficients, say, integer, floating-point or imaginary numbers, it's also important that expression trees, also implemented as objects, can easily be used as symbols, for example, of the polynomial object. But it is not known in advance what the expression tree of a particular system will look like.

The Term Object

These considerations lead to defining a user-level term object, as follows :

@interface Term : Object
{
	id		coefficient;
	id		symbol;
	int		exponent;
}

The OBJECTIVE-C class declaration above defines coefficients and symbols to be arbitrary OBJECTIVE-C objects (the C-type id ), while we choose for exponents that are small C integers (the C-type int ).

If the exponent is equal to zero, then the symbol of the term is the nil object (no object). As part of a polynomial, it definitely makes sense for a term to have an exponent equal to zero, while having a non-nil symbol. But for terms as stand-alone objects, it is better to have a unique representation of the constant term.

The decision to make exponents small (and, to be used with the Polynomial object, non-negative) integers, instead of using generalized object exponents, was based, partly on the fact that object exponents are usually dealt with at the Expression level, and partly on the conviction that choosing for the more general solution would go in this case at the expense of the intuitivity and simplicity of the Polynomial object : it would have made giving a definition of in-exact polynomial division, for instance, more difficult. And it is known how to use simple polynomials as building-blocks for constructing Laurent or Puisseux polynomials (with negative or fractional exponents respectively). We felt that it would be safer to provide separate classes for these mathematical objects, with sensible implementations of the associated elementary arithmetical operations.

Another ``strategical'' decision was, whether to define a Term as a triple of coefficient, symbol and exponent, or whether to define a Term as a pair of coefficient and auxiliary Power object. The latter object was in earlier versions of the Computer Algebra Kit then defined as a pair of symbol and exponent. We finally decided against this; Power objects were not really useful. Simple operations like adding or differentiating already make it necessary to convert towards the notion of Term . And a Power object can be considered as a Term with a coefficient equal to one (a monic term).

By adopting the current definition of Term , it is possible to define recursive polynomials by using just two abstractions -- instead of three -- which is a significant simplification of the surface area of the library i.e., the number of notions that a library user has to get familiar with, and that should not be underestimated. Whenever possible, we try to reduce the number of objects that the application programmer has to get acquainted with.

Recursive Polynomials

Recursive polynomials -- we are mostly following the terminology of Zippel's book [Zippel, 1993] -- are abstractly sequences (sums) of Term objects. The coefficients of the terms of such polynomials, can be recursive polynomials again, in a variable (symbol) less i.e., the recursion is on the number of variables.

For example, here's a polynomial in two variables written in a recursive form :

(x2 + 1) y2 + x + 1

The coefficient of the term in y2 is a polynomial in x.; not all coefficients are simply integers.

However, requiring the user of the Polynomial object to fix the number of variables, is not, in general, a good way to define when the recursion ends. It would have been a plausible choice for variable dense recursive polynomials, but not for variable sparse recursive polynomials. The former kind of polynomials are exactly defined by a base scalar domain (integers, in the example above) and a set of admissible symbols (x and y for example). If some variable is ``missing'', then a power with exponent zero is substituted. For example, in the variable dense representation, the polynomial above has two terms :

(x2 + 1) y2 + (x + 1) y0

Variable sparse polynomials on the other hand dynamically adjust their set of symbols. Only symbols with a non-zero exponents count. The example polynomial has three terms, as a variable sparse recursive polynomial.

The dense variety of recursive polynomials is useful for a class of algorithms, including the subresultant and polynomial greatest common divisor ones. The sparse kind is much more convenient when doing substitutions etc., because one can freely mix polynomials in all sorts of symbols.

Either way, for recursive polynomials in the Computer Algebra Kit, the kernel implementor has to fix the base scalar domain of the Polynomial object being used. The innermost coefficient objects of a recursive polynomial are scalar objects. Perhaps it would have been possible to simply consider any non-polynomial object as scalar, because in general scalar objects will be Integer , Fraction , Complex , Float or perhaps even some kind of Expression object, but that would have excluded the possibility of working with multivariate polynomials over a polynomial base scalar domain, which also might be interesting in some cases.

In other words, a Polynomial is created using the newScalar: method, which takes an object, such as an Integer object, as argument and then returns a recursive, variable sparse polynomial object with integer base scalar coefficients :

@interface Polynomial : Object { ... }

+ newScalar:aScalar;

- makeVariableDense;
- makeVariableSparse;

- (BOOL)isVariableDense;
- (BOOL)isVariableSparse;

- eachTerm;
- insertTerm:aTerm;
- removeTerm;

@end

It is then possible to use any of the insertTerm: and removeTerm methods to insert and remove terms, or to obtain a sequence of the terms in the polynomial using eachTerm. Like this, a variable sparse polynomial in any number of symbols can be constructed.

A trick that is very useful to reduce the number of public methods in the interface of the Polynomial object, is to use a given representation as starting-point for creating another one.

Given the variable sparse polynomial, the user can obtain the equivalent polynomial in the dense representation using the makeVariableDense method. This method will collect all symbols that occur in the sparse polynomials, and use that particular set of symbols for creating the corresponding variable dense polynomial.

This is certainly not equivalent to enumerating all sorts of constructors in the public interface for all of the possible representations. Methods with many arguments are confusing; mechanisms for inserting ``default'' values in the constructors only make things worse. And listing many methods with similar names and only slight differences in functionality make it difficult for a user (and for the implementor after a few years !) to decide which one to use.

The goal when designing the public interface of the Polynomial object, was to choose short, descriptive method names, and a means for doing so, was to use polynomial representations that are easier to construct, as intermediate steps for creating the ones that are easier to compute with.

A related issue is, whether there should be just one Polynomial object in the library, or a whole range of classes, for different representations. The continuous conversion from one representation to another, to have the ``right'' representation for the problem at hand, and the need for making different representations smoothly complement each other -- the notion of scalar for a variable sparse polynomial must be identical to that of a scalar for a variable dense polynomial ! -- , definitely justifies having just one Polynomial object, and considering polynomial representation as a run-time, instead of a compile-time attribute of the polynomial. Besides, using different classes for all the representations would have dramatically increased the size of the interface of the library.

Degree Dense and Degree Sparse Polynomials

As said before, the library user decides whether terms with exponent equal to zero occur in the sequence of terms, or not. A similar choice must be made for the coefficients : in the degree sparse representation of a polynomial, only terms with a non-zero coefficient are significant. In the degree dense representation the ``gaps'' between non-adjacent exponents are filled with terms with a coefficient equal to zero.

For recursive polynomials, it is easy to order the terms in the sequence associated to the polynomial. The leading term is the term with the greatest exponent; it is the first term in the sequence returned by eachTerm and it always has a non-zero coefficient. The degree of a recursive polynomial is equal to the exponent of the leading term. If the polynomial is equal to zero, then we set the degree to minus one, since we consider a zero polynomial to be identical to an ``empty'' polynomial i.e., a polynomial containing no terms. A non-zero constant polynomial has degree equal to zero. It is exactly for this kind of thing that we prefer typing quantities like degrees and counts to the C type int instead of unsigned. The order of the polynomial is set to the exponent of the last proper term.

@interface Polynomial : Object { ... }

- makeDegreeDense;
- makeDegreeSparse;

- (BOOL)isDegreeDense;
- (BOOL)isDegreeSparse;

- (int)order;
- (int)degree;
- (int)numTerms;

@end

By default, newScalar: creates a degree sparse, recursive and variable sparse polynomial. It is not necessary to provide any other constructors, because the method makeDegreeDense returns an equivalent degree dense polynomial. The sequence returned by eachTerm will now contain padding zero terms.

Degree dense polynomials are convenient to use in the implementation of, for example, univariate polynomial factorization algorithms.

The Monomial Object

A monomial consists of a scalar object multiplied by a sequence (product) of monic terms i.e., powers of symbols. For example, the following is a monomial :

3 x2 y4 z2

Again, both scalar and symbols of a monomial can be arbitrary objects.

Terminology in this area is not really settled. Some people don't include the scalar coefficient in the definition of a monomial, and call term what we call monomial. Our choice is motivated by the fact that we just need a single concept, Term , -- that we can even share with the recursive polynomials -- for defining monomials, and just two concepts, Term and Monomial , for defining expanded polynomials.

The definition of a Monomial parallels to some degree that of a recursive polynomial, except for the fact that we're working now with a product instead of a sum of terms, and because of an additional scalar factor :

@interface Monomial : Object
{
	id		scalar;
	id		terms;
}

+ newScalar:aScalar;

- makeVariableSparse;
- makeVariableDense;

- scalar;
- eachTerm;
- insertTerm:aTerm;
- removeTerm;

- (int)compareTerms:aMonomial;

@end

The interface doesn't export the collection of terms of a Monomial ; this collection can only be accessed by means of the eachTerm method, as a Sequence. Doing otherwise would have forced us to define an additional PowerProduct object, which we want to avoid at all events, to reduce the size of the library interface.

Also note that, just like in the case of recursive polynomials, there is a choice for monomials on whether to fix or not the set of allowable symbols, and whether or not to store symbols with exponent equal to zero. The issues mentioned before with respect to constructors, conversions etc. apply here as well.

There exists a whole mathematical theory on ordering monomials, comparing them to each other. We will not discuss that here; it suffices to say that, a run-time attribute of the Monomial object is the way it is ordered, and that by default we use the lexicographic ordering. For that ordering we have, if x is smaller than y, and if y is smaller than z :

x2 y4 z2 < x2 y5 z

For variable dense monomials, the notion of lexicographic ordering is defined by the set of symbols of the monomial, while for a variable sparse monomial, created by newScalar:, it's the implementation of compare: of the Symbol class that counts.

Expanded Polynomials

Expanded polynomials are abstractly sums of Monomial objects. The monomials are ordered with respect to some ordering. The sequence returned by eachMonomial has, as its first member, the monomial that is largest with respect to that ordering. The degree of an expanded polynomial is the degree (sum of exponents) of this leading monomial.

@interface Polynomial : Object { ... }

- makeExpanded;
- makeRecursive;

- (BOOL)isExpanded;
- (BOOL)isRecursive;

- eachMonomial;
- insertMonomial:aMonomial;
- removeMonomial;

@end

The notion of base scalar coincides in the recursive and expanded case. For example, the recursive polynomial :

(x2 + 1) y2 + x + 1

can be written as an expanded polynomial like :

x2 y2 + y2 + x + 1

In both cases the base scalar ring consists of Integer objects. Expanded polynomials are created by using the makeExpanded method. Just as in the recursive case, there are some variations possible, such as, whether the expanded polynomial is degree dense or degree sparse, and whether it is variable sparse or variable dense.

Expanded polynomials are commonly used in polynomial standard basis algorithms.

The Files Polynomial.h and Polynomial.rtf

The programmer that implements a kernel of a symbolic computation system, knows what kind of functionality the Polynomial building block features, by consulting the header file of the object. This file should also be included in the source file, when referring to the object somewhere in an OBJECTIVE-C program. Many of the code snippets given so far, were extracted from this header file Polynomial.h.

The Polynomial object implements more operations than listed in the public header file; these undocumented methods are private. Working with private methods and private objects is essential for hiding implementation details to the consumer of the library.

In the spirit of the book of Kernighan and Pike [Kernighan and Pike, 1984], that advocates building small tools for solving this sort of work, we wrote a UNIX filter (also implemented in OBJECTIVE-C) that makes it possible to list some OBJECTIVE-C methods as public and others as private . During compilation of the Computer Algebra Kit, one works with the private header file, while, once the object is compiled, only the public header file counts. Both header files, and the source file of the object, are extracted from a single source file.

Another option of the same UNIX filter, extracts an OBJECTIVE-C specification sheet from the source file. Such a specification sheet doesn't simply list the methods of the Polynomial object, but it also documents the return value, side effects, and arguments of each method, and how the object as a whole fits in the framework of the Computer Algebra Kit.

Private Polynomial Objects

There are currently, in addition to the Polynomial object, eight special purpose objects for different polynomial representations. These objects carry names such as,

@interface VarDenseRecursiveDegDensePolynomial : Object {}
and they are -- for obvious reasons ! -- private.

For these objects, the polynomial representation is fixed at compile-time. They are meant to implement the low-level mapping of abstract methods, such as insertTerm:, towards a concrete data-structure, such as an array of coefficients, or a linked list of C terms (not to be confused with the public OBJECTIVE-C Term object).

The Polynomial object then, implements some methods, by making the carrier object -- the private object corresponding to the current internal representation -- do the real work :

- insertTerm:aTerm
{
	if ([self isRecursive]) {
		[carrier insertTerm:aTerm];return self;
	} else {
		return [self error:"Not recursive"];
	}
}
This class of private objects typically only performs simple, linear, O(n)-style operations : inserting a term or monomial into the list or array of terms or monomials, adding two polynomials together, multiplying a polynomial by a scalar or term etc.

Multiplication, division, factoring, gcd-computation of polynomials etc. is implemented at the Polynomial object level. This follows of course from the observation that, in polynomial algorithms, the bulk of the time is spent in doing very low-level operations, like making linear combinations of polynomials.

Producing Private Objects

The eight different polynomial classes for specific polynomial representations are further specialized into a range of subclasses that correspond to different base scalar or coefficient domains. For example, objects such as :

DoubleVarDenseRecursiveDegDensePolynomial
or,
ObjectVarDenseRecursiveDegDensePolynomial
are two distinct, internal OBJECTIVE-C objects for working with polynomials with floating-point (of C type, double ) or OBJECTIVE-C objects (of C type, id) in the particular representation. For these classes, both scalar type and polynomial representation is fixed at compile time.

The polynomial object with double coefficients is implemented as :

struct dbl_vardnsrecdegdns_poly { 
	int		numTerms;
	double *	coefficients;
}
while the polynomial object with id coefficients is essentially implemented as :

struct obj_vardnsrecdegdns_poly {
	int		numTerms;
	id *		coefficients;
}
In other words, when working with a variable dense, recursive and degree dense Polynomial object with Float objects as base scalars, the system selects the corresponding internal representation, and when, for example, adding such polynomials together using the method add:, then the work is done , in the end, by plain statically, tightly bound C code :

dbl_vardnsrecdegdns_add(c,a,b,n)
	double *c,*a,*b;
	int n;
{
	while (n--) { *c++ = (*a++) + (*b++); }
}
Currently there are such internal, special purpose polynomial objects for polynomials with integer, modular and floating point coefficients, and this can be extended in the future, without changing the public interface of the Polynomial object. When working with other kinds of coefficients than the one listed above, the system can fall back on the representation for the id type, which stands for arbitrary, loosely bound, OBJECTIVE-C objects :

obj_vardnsrecdegdns_add(c,a,b,n)
	id *c,*a,*b;
	int n;
{
	while (n--) { *c++ = [(*a++) add:(*b++)]; }
}
The Polynomial class can therefore work over coefficient domains that are defined by the kernel implementor, and that are not part of the Computer Algebra Kit.

On the other hand, for many polynomial coefficient domains -- rational function fields for example -- it simply makes no sense to replace the object representation by a special purpose one, and to tightly bind the coefficient domain to the polynomial algorithms : the overhead of Smalltalk-like run-time binding is in these cases neglectible when compared to the cost of arithmetic in the particular coefficient domain.

Following the over-all philosophy of [Kernighan and Pike, 1984] -- designing auxiliary development tools to manage the software -- we implemented yet another precompiler for producing the dozens of low-level, repetitive and simple-minded OBJECTIVE-C classes that are the building bricks for implementing the Polynomial object. This tool provides a statement similar to the C preprocessor statement :

include "package.h" (params)
, but it accepts parameters in the style of ADA-style packages for file inclusion with parametrized types. The tool is prepended to the compiler chain that is used for compiling the Computer Algebra Kit.

This is sufficient for our purposes; because we limit ourselves to applying static-binding techniques on the O(n) algorithms. At that level, there is no need for incremental compilation, and no memory allocation difficulties, that are typical when trying to extend early-binding to higher-level algorithms. Incremental compilation and uniform memory allocation (all OBJECTIVE-C objects are allocated in a ``heap'') are benefits of the late-binding approach.

Conclusion

Put in another way, parametrized packages and plain vanilla, statically bound C code are used in the Computer Algebra Kit as a production technology at the lowest level for a few time-critical loops, well hidden from the user, while objects, Smalltalk-like objects, are the consumer-oriented technology and the top-level interface that we want to present to the kernel programmer.

From a mathematical point of view, the Polynomial object allows one to compute over more coefficient domains, over dynamically constructed domains, in more representations etc., compared to corresponding function package technology. For example, the square free factorization algorithm of the Computer Algebra Kit, decomposed polynomials over the Gaussian integers and over Galois fields, or handled polynomials with coefficients in a function field. An implementation of a standard basis algorithm was able to compute with polynomials with coefficients in an integral domain, such as the integers, the Gaussian integers etc. or in a field, such as the rational numbers, the modular integers etc. An implementation of the Sturm algorithm, for counting real roots of a polynomial, works for degree dense and degree sparse polynomials the like, over ordered domains. The Cantor-Zassenhaus algorithm factored polynomials into irreducible parts over the integers modulo a small prime and over Galois fields.

The Polynomial object is, more importantly perhaps, from the software engineering or system-building point of view, some intermediate level between the top-level of the kernel of a symbolic computation system, and the low-level function packages that in general form the implementation of a kernel. Grouping together a range of function packages, it is the key component in the implementation of the Expression object of a small, extensible interactive system that comes with the Computer Algebra Kit. As a reusable, loosely coupled component, it is compatible with objects from other object libraries, such as objects for transferring objects over a network or such as objects for storing objects on disk. It is also compatible with, and has been used with, a range of exisiting development tools, including optimizing and incremental OBJECTIVE-C compilers, browsers, debuggers and profilers, that are of great value when implementing or evaluating some particular algorithm.

References

Cox, 1986
Cox, B. J. (1986).
Object-Oriented Programming : An Evolutionary Approach.
Addison-Wesley, Reading, Mass.

Goldberg and Robson, 1983
Goldberg, A. and Robson, D. (1983).
Smalltalk-80 : The Language and its Implementation.
Addison-Wesley, Reading, Mass.

Kernighan and Pike, 1984
Kernighan, B. W. and Pike, R. (1984).
The Unix Programming Environment.
Prentice-Hall, Inc.

Zippel, 1993
Zippel, R. (1993).
Effective Polynomial Computation.
Kluwer Academic Publishers, Boston, Dordrecht, London.



Footnotes

...We
David Stes, Molenstraat 5, 2018 Antwerpen, Flanders, E-mail: stes@pandora.be

...years
The work was done independently, and we would like to thank our family and friends for the support over the years.



David Stes
6/10/1998