Skip to content


My polynomial library provides a Python class for manipulating one variable polynomials with rational coefficients.   It was originally based on the ActiveState polynomial library published by Sam Denton but is now significantly different.  More recently I have taken on board ideas exhibited in Jim Randell’s own polynomial class in his library.

The library was originally first used in the Funny Dice teaser published by Michael Fletcher in the Sunday Times.   Mike has recently provided another teaser on the same theme with Wonky Dice which moves from pairs of six sided dice to octagonal ones.  Here is a general routine for solving dice probability puzzles of this type:

This depends on code to generate cyclotomic polynomials which I first implemented as follows:

This is fairly fast for low \(n\) values but runs out of speed for \(n\) in the hundreds. Jim Randell’s version here is slower for low \(n\) but has the advantage that it can be cached and this makes it much faster for large \(n\).

In looking at various faster approaches, I came across a brilliant piece of Python code for producing cyclotomic polynomials developed by David Furnish (this is on the Rosetta Code site). David dispenses with conventional polynomial arithmetic and expresses cyclotomic polynomials in the form: \[\large \Phi_n(x)=\frac{x^n-1}{\prod_\limits{d|n\;d\lt n}\Phi_d(x)}\]

This can be applied recursively, which means that all evaluations will eventually reach cyclotomic polynomials for \(d\) equal to a prime, which can be easily generated: \[\large\Phi_p(x)=\sum_{n=1}^{p-1}x^n\;=\;\frac{x^p-1}{x- 1}\]From the above formulae we can see that all the polynomials involved in computing cyclotomic polynomials are of the form \((x^p – 1)\), which allows a series of operations on these polynomials to be represented as list of integers where an integer \(n\) represents the polynomial \((x^n – 1)\). This allows a cyclotomic polynomial to be represented in a very compact notation, for example:\[C(305)=\frac{(385,5,7,11)}{(35,55,77,1)}\] If we have a polynomial expressed as an ordered sequence of (coefficient, power) pairs and we want to multiply it by \((x^p – 1)\) (represented by the integer \(p\)), this will produce two sequences of pairs, one with its coefficients negated and a second with its powers increased by \(p\). These two streams can then be merged to give the result. On the other hand if we wish to divide by the term instead of multiplying, we can convert the input (coefficient, power) stream into one that has the powers reduced by \(p\) whilst also merging this term back into the input stream so that it will cancel this term’s negated twin when the merge is performed. This all works because we know that the division is an exact one with no remainder.

David Furnish implements this approach using Python generators to produce a blindingly fast way of producing cyclotomic polynomials of large degree. To illustrate its speed we can look for the first cyclotomic polynomials that have coefficients of a particular integer magnitude. For example, C(105) is the first that contains a coefficient of magnitude 2; the first coefficient of magnitude 3 appears in C(385).

This program times finding these results using the David’s approach which I now use in my polynomial library:

and here are the results: