Polynomials

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 enigma.py 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:


from itertools import product
from polynomial import Poly, cyclotomic

# find dice combinations giving the same probabilties as
# standard dice pairs
def sicherman(n):

  # return the dice represented by a polynomial
  def die(ct):
    if sum(ct) == n and all(c <= 0 for c in ct):
      return sum(((i,) * c for i, c in enumerate(ct, 1)), ())

  # form the polynomial representation of a fair octagonal
  # die as the product of cyclotomic polynomials (less the
  # x factor)
  # (see https://en.wikipedia.org/wiki/Sicherman_dice)
  fctrs = []
  for d in range(2, n + 1):
    if n % d == 0:
      fctrs.append(cyclotomic(d))
  
  # consider all subsets of the factors of this polynomial
  for ms in product((0, 1, 2), repeat=len(fctrs)):
    p1, p2 = Poly([1]), Poly([1]) 
    for m, f in zip(ms, fctrs):
      mp = (1, f, f * f)
      p1 *= mp[m]
      p2 *= mp[2 - m]
    d1 = die(p1)
    if d1:
      d2 = die(p2)
      if d2 and d1[-1] <= d2[-1]:
        yield d1, d2

from itertools import product

from polynomial import Poly, cyclotomic

# find dice combinations giving the same probabilties as

# standard dice pairs

def sicherman(n):

# return the dice represented by a polynomial

def die(ct):

if sum(ct) == n and all(c <= 0 for c in ct):

return sum(((i,) * c for i, c in enumerate(ct, 1)), ())

# form the polynomial representation of a fair octagonal

# die as the product of cyclotomic polynomials (less the

# x factor)

# (see https://en.wikipedia.org/wiki/Sicherman_dice)

fctrs = []

for d in range(2, n + 1):

if n % d == 0:

fctrs.append(cyclotomic(d))

# consider all subsets of the factors of this polynomial

for ms in product((0, 1, 2), repeat=len(fctrs)):

p1, p2 = Poly([1]), Poly([1])

for m, f in zip(ms, fctrs):

mp = (1, f, f * f)

p1 *= mp[m]

p2 *= mp[2 - m]

d1 = die(p1)

if d1:

d2 = die(p2)

if d2 and d1[-1] <= d2[-1]:

yield d1, d2

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


from functools import reduce
from operator import mul
from number_theory import factors
from polynomial import Poly

def cyclotomic(n):
  '''
  Return the n'th Cyclotomic Polynomial

  >>> print(cyclotomic(7))
  x^6 + x^5 + x^4 + x^3 + x^2 + x + 1

  >>> print(cyclotomic(15))
  x^8 - x^7 + x^5 - x^4 + x^3 - x + 1

  >>> print(cyclotomic(30))
  x^8 + x^7 - x^5 - x^4 - x^3 + x + 1

  >>> t = cyclotomic(105); t[:10]
  (1, 1, 1, 0, 0, -1, -1, -2, -1, -1)
  '''
  # phi(1) = x - 1
  phi = Poly([-1, 1])
  if n == 1:
    return phi
  # find the distinct prime factors of n (i.e. f is square-free)
  f = set(factors(n))
  
  # If p is a prime that does not divide m then:
  #
  #    phi(m.p)[x] = phi(m)[x^p] / phi(m)[x]
  #
  for p in f:
    phi, rem = divmod(phi.compose(Poly.from_cepairs([(1, p)])), phi)
  
  # if q is the product of the square-free factors of n:
  #
  #    phi(n)[x] = phi(q)[x^(n/q)] 
  #
  return phi.compose(Poly.from_cepairs([(1, n // reduce(mul, f, 1))]))

from functools import reduce

from operator import mul

from number_theory import factors

from polynomial import Poly

def cyclotomic(n):

'''

Return the n'th Cyclotomic Polynomial

>>> print(cyclotomic(7))

x^6 + x^5 + x^4 + x^3 + x^2 + x + 1

>>> print(cyclotomic(15))

x^8 - x^7 + x^5 - x^4 + x^3 - x + 1

>>> print(cyclotomic(30))

x^8 + x^7 - x^5 - x^4 - x^3 + x + 1

>>> t = cyclotomic(105); t[:10]

(1, 1, 1, 0, 0, -1, -1, -2, -1, -1)

'''

# phi(1) = x - 1

phi = Poly([-1, 1])

if n == 1:

return phi

# find the distinct prime factors of n (i.e. f is square-free)

f = set(factors(n))

# If p is a prime that does not divide m then:

# phi(m.p)[x] = phi(m)[x^p] / phi(m)[x]

for p in f:

phi, rem = divmod(phi.compose(Poly.from_cepairs([(1, p)])), phi)

# if q is the product of the square-free factors of n:

# phi(n)[x] = phi(q)[x^(n/q)]

return phi.compose(Poly.from_cepairs([(1, n // reduce(mul, f, 1))]))

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:


  # find first occurence of coefficient 'cf'
  from time import process_time_ns as ns_cntr
  t0 = ns_cntr()
  n, cf = 0, 1
  while cf < 10:
    n += 1
    cp = cyclotomic(n)
    if cf in cp or -cf in cp:
      ix = cp.index(cf) if cf in cp else cp.index(-cf) 
      t = ns_cntr() - t0
      print(f'C[{cf}]: {n:4} @{ix:6} {t / 1000000:8.3f}ms')
      cf += 1

# find first occurence of coefficient 'cf'

from time import process_time_ns as ns_cntr

t0 = ns_cntr()

n, cf = 0, 1

while cf < 10:

n += 1

cp = cyclotomic(n)

if cf in cp or -cf in cp:

ix = cp.index(cf) if cf in cp else cp.index(-cf)

t = ns_cntr() - t0

print(f'C[{cf}]: {n:4} @{ix:6} {t / 1000000:8.3f}ms')

cf += 1

and here are the results:


C[1]:    1 @     1    0.000ms
C[2]:  105 @     7    0.000ms
C[3]:  385 @   119   31.250ms
C[4]: 1365 @   196  265.625ms
C[5]: 1785 @   137  437.500ms
C[6]: 2805 @   588 1031.250ms
C[7]: 3135 @   616 1281.250ms
C[8]: 6545 @  1875 5250.000ms
C[9]: 7917 @  1753 7656.250ms

C[1]: 1 @ 1 0.000ms

C[2]: 105 @ 7 0.000ms

C[3]: 385 @ 119 31.250ms

C[4]: 1365 @ 196 265.625ms

C[5]: 1785 @ 137 437.500ms

C[6]: 2805 @ 588 1031.250ms

C[7]: 3135 @ 616 1281.250ms

C[8]: 6545 @ 1875 5250.000ms

C[9]: 7917 @ 1753 7656.250ms

Polynomials

Recent Posts

Recent Comments

Archives

Categories

Meta