In the late 80s/early 90s I had developped a simple Symbolic algebra in Smalltalk-80 (v2.3-v2.5 from Parcplace).
One of the major challenges was to implement algorithms for univariate and multivariate polynomials, over a commutative Ring (mostly the rationals).
I wanted to redo this work in Squeak Smalltalk for a long time.
So here we go, there is now a first version at squeaksource:
http://www.squeaksource.com/STEM/Polynomials-nice.1.diff
(well that was the early version when I started writing this post, take the head revision for better experience)
There are a few prerequisites that can be found in STEM like IntegerFactorization and CollectionArithmetic
Some design choices deserve comments.
First let's review the vocabulary.
A univariate polynomial $p$ is like a function of one variable, let's note it
$p(x_1)$, involving a sum of powers of the variable, up to some degree $n$, multiplied by some coefficients $c_i$,
which are elements of a ring $R$.
\[ p(x_1) = \sum_{i=1}^{n} c_i x_1^i \: , \: c_i \in R \]
For example, with integer coefficients $c_i \in \mathbb{Z}$:
\[ p(x_1)=1 + 2 x_1^3 \]
A multivariate polynomial is like a function of several variables, say $x_1, \ldots , x_r$
For example
\[ p(x_1,x_2,x_3)=1 + 3 x_1 x_2 - 5 x_1^2 x_3 + 2 x_1 x_2^2 x_3 - x_3^5\]
We can see a polynomial as a sum of terms.
Each term is the product of a coefficient in underlying Ring, $c \in R$ and a monomial.
We call it a monomial coefficient.
A monomial is the product of the variables, each raised to some powers.
Here we have those terms:
$$+1 \: x_1^0 x_2^0 x_3^0$$
$$+3 \: x_1^1 x_2^1 x_3^0$$
$$-5 \: x_1^2 x_2^0 x_3^1$$
$$+2 \: x_1^1 x_2^2 x_3^1$$
$$-1 \: x_1^0 x_2^0 x_3^5$$
Each monomial is represented by anArray, where (anArray at: i) represents the power of variable of rank i.
Somehow, we can represent the polynomial with an Array of monomial coefficients, and the associated monomial powers.
For the example above, that would be:
- monomial coefficients:
#(1 3 -5 2 -1) - monomial powers:
#( #(0 0 0) #(1 1 0) #(2 0 1) #(1 2 1) #(0 0 5) )
That's how multivariate polynomials are represented in above package.
Note that the variables rank spans 1 to: 3, the size of each monomial power.
If there is only, say powers of $x_2$ and $x_{99}$, and the degree vs any other variable is zero, then that makes this storage strategy very inefficient.
That's why a SparseArray is used instead, above a certain rank.
See also NotSoSparseArray which is just an Array understanding the SparseArray specific protocol.
For univariate polynomials, the monomials involve the power of a single variable.
Thus $p(x_1)=1 + 2 x_1^3$ is represented with:
- monomial coefficients:
#(1 2) - the associated powers:
#(0 3) - and variable rank:
1
This variable rank is required because it cannot be encoded in the powers unlike MultivariatePolynomial
and because we want to be able to produce a multivariate polynomial by operating on univariate ones like $p(x_1,x_2)=q(x_1) \cdot r(x_2)$
Another design choice concerns the polynomial of degree zero. We could just use an element of underlying ring $R$, that would be mostly an instance of Number. But if we do this while still wanting to perform polynomial requests like telling about its degree, leading coefficient, or any other polynomial stuff, that would mean implementing the whole polynomial protocol in Number.
That sounds too much, the Number methodDictionary is already bloated.
Instead, we design a special PolynomialDegree0 and must use 2 asPolynomial if we want to consider 2 as a polynomial.
And we also design Polynomial classes that perform automatic reduction, like a Fraction with unit denominator reduces to an Integer:
- a
MultivariatePolynomialdoes automatically reduce to aUnivariatePolynomialwhen only one variable has degree greater than 0 - Both
MultivariatePolynomialandUnivariatePolynomialdo reduce automatically toPolynomialDegree0when their total degree is zero
Next design choice concerns the ordering of monomials:
- a
UnivariatePolynomialalways fuse monomials of same power - a
UnivariatePolynomialalways order its monomials in descending powers - a
MultivariatePolynomialalso fuse monomials of same powers in each variable - a
MultivariatePolynomialalways order its monomials in descending powers of last variable, then last but one, etc...
All the invariants that should be enforced can be found in class comments.
With that rules, we have a canonical representation of polynomials: every two equivalent polynomials will share the same representation.
That's a nice feature for testing polynomial equality.
But we soon break this rule: one of the major operations we want to be able to perform on polynomials is to factor them into a product of irreducible terms.
This reduction is not cheap, and we would like to preserve the factorized structure, so we soon need more classes in the Polynomial zoo:
- a
PolynomialPowerfor raising a polynomial to a positive integer power - a
PolynomialProductfor representing the product of several polynomials, or at least a polynomial and a scalar
There are rules for automatically reducing these instances when they degenerate, like raising to power 0 or 1, and multiplying a list of polynomials of size 0 or 1, and for some of the trivial combinations like:
- a power of power $(p(x)^m)ˆn$ is reduced to a flat power $p(x)^{m \cdot n}$
- a power of product $(p(x) \cdot q(x))^n$ is reduced to a product of powers $p(x)^n \cdot q(x)^n$
- a product of product $(p(x) \cdot q(x)) \cdot r(x)^n$ are flatten $p(x) \cdot q(x) \cdot r(x)^n$
- a product of powers $p(x)^m \cdot p(x)^n$ is reduced to a flat power $p(x)^{m + n}$
But it is no more a canonical representation, nothing prevents using a product with partial factorization, that is where some terms would still be reducible.
One way to come back to a more canonical form is to either:
- completely factor the polynomial with
p factor
though it does not give any guaranty about the order of terms - completely expand the polynomial with
p expandPolynomial
To check if two different representations p and q are the same polynomial, we can use the query p isSamePolynomialAs: q.
Note however that the implementation is not particularly efficient.
So what can we do with those polynomials?
First create the variables:
Then of course, performing arithmetic like any other algebraic object in Smalltalk.
Note below how the polynomials can perform mixed arithmetic with numbers, no need to use asPolynomial at this stage:
With Smalltalk precedence this is just the Horner form $((x_1+2) \cdot x_1-1) \cdot x_1^2+4$
If preferring the monomial form at the price of parenthesis here is another one:
Then we can inquire about degree, coefficients, etc..
Since a polynomial can be viewed as a function, we can derivate and get $\frac{\partial t(x_1,x_2)}{\partial x_1}$
Or just evaluate it:
For composing, $p \circ q$, one can simply use evaluation:
,
or a bit more efficient:
Beware, however, evaluation of MultivariatePolynonial currently only works for numerical argument.
Other arithmetic we can perform is Eclidean division, greatest common divisor or least common multiple.
Or we can inquire a factored form:
There is more we can do like p discriminant, t resultant: p
And finally we can also play with rational functions (division of polynomials) which more or less work like Squeak Fraction
with a rather higher cost of course since it involves gcd and lcm computations.
There are limitations though. The implementation is only effective for the ring of integers $\mathbb{Z}$ and the field of rationals $\mathbb{Q}$.
It would only works on commutative ring, anyway, for example not on a ring of square Matrices nor on Quaternions
The implementation will also very hardly work with floating point, because floating point do not have the properties of a ring with + and *.
Floating point operations are indeed not associative to begin with.
For example, evaluating a polynomial with Float coefficients thru p value: 0.1 or p valueHorner: 0.1 would probably not even lead to the same value.
Since the multivariate polynomials $R[x_1,\dots,x_n]$ over the ring $R$ is also a ring, it can be viewed as an univariate polynomial in $x_n$ over the ring $R[x_1,\dots,x_{n-1}]$, that is $R[x_1,\dots,x_{n-1}][x_n]$ or recursively, $R[x_1]\dots[x_n]$
My historical implementation in Smalltalk-80 did use such a representation where coefficient of polynomial in $x_n$ were polynomials in $x_{n-1}$ and recursively.
Current implementation would probably not work either if using polynomials as monomialCoefficients, the arithmetic operations would at best bring the $R[x_1,\dots,x_n]$ representation back.
If you are interested in a more general symbolic algebra making the algebraic structures explicit, generalizing their properties and generalizing the algorithms, then the notion of domain becomes essential, like explored by IBM scratchpad (axiom) in the 70s.
In this case, I must strongly recommand to inquire the remarkable work of Luciano Notarfrancesco at https://github.com/len/Arrows
Final notice: the code has not been generated by AI. Nor this post.
However, while writing this post, I did massively question Gemini, either for dumb things like remembering some LaTeX snippets that I did not use for 20 years, or for verifying some mathematical properties, because it's so much more efficient than traditional internet search.
While doing so, this gave me many new ideas I'd like to explore, if life isn't too short.
I'm pretty sure that I could ask an AI to write this entire blog post, I am not mentally ready to do this now, since I don't see the point, and my prose is already quite verbose, no need to worsen that style.
