"""Polynomials. Coefficients are assumed numeric. Only natural powers are considered. $Id: polynomial.py,v 1.19 2007/06/03 16:42:30 eddy Exp $ """ import types from study.snake.lazy import Lazy class invalidCoefficient (TypeError): "Invalid coefficient for polynomial" class unNaturalPower (TypeError): "Power of variable in polynomial is not a natural number" class Polynomial (Lazy): """Model of the mathematical ring of polynomials (in one free variable). Supports arithmetic (+, -, *, /, %, divmod, power, negation), comparison (but cmp() will yield zero in some cases where == will say no), representation (as a python lambda expression), evaluation (i.e. calling as a function), repeated integration (as <<) and differentiation (as >>), use as boolean (only the zero polynomial is false) and (lazy) hashing. Lazy attributes: =============== rank -- highest power of the free variable with non-zero coefficient normalised -- same polynomial scaled so rank's coefficient is 1 derivative -- result of differentiating the polynomial assquares -- decompose as sum of scaled squares and a remainder sign -- like cmp(self, 0) but None when ill-defined or variable isreal -- are all coefficients real ? factors -- tuple of normalised irreducible factors and an optional scalar Note that a poly with true isreal considers positive-definite quadratic factors to be irreducible; assigning isreal = None will persuade a poly to believe quadratic factors are reducible. The value of p.assquares is of form (bok, rem) with rem either zero or of odd degree and p-rem equal to reduce(lambda y, (x, v): y + x*x*v, bok.items(), 0). The value of sign may presently sometimes be None when it needn't be; I've yet to find an example, but the computation in use only deals with `easy enough' cases. Note that cmp() yields zero when the difference's sign is None or zero; but == is true only when sign is zero (actual equality). Methods: ======== coefficient(n) -- yields coefficient of (: x**n ←x :) in polynomial hcf([poly, ...]) -- yields highest common factor of arbitrarily many integral([start=0, base=0]) -- integration unafter(poly) -- input to poly yielding self as output seek_root([guess=0, tol=1e-6]) -- find an input mapped to zero seek_factor([guess=None]) -- find a factor, if possible See individual methods' docs for details.\n""" def __init__(self, *args): """Constructor. If exactly one argument is passed and it is a mapping, its keys are used as powers and values as coefficients of each power's term in the polynomial; if exactly one argument is passed and it is a sequence, it is interpreted as a mapping whose keys are indices into the list and values are entries in the list. Otherwise, each argument must be numeric and the arguments are interpreted as if they'd been passed as a sequence. Thus: Polynomial({ power: coeff ... }) == (: sum(: coeff * x**power :) <- x :) Polynomial([ co_0_, co_1_ ... ]) == Polynomial(( co_0_, co_1_ ... )) == Polynomial( co_0_, co_1_ ... ) == (: sum(: co_i_ * x**i <- i :) <- x :) As special cases of the last, Polynomial( value ) generates the constant polynomial value <- x and Polynomial() generates the zero polynomial - but note that this is the scalar zero polynomial; if a vector zero or dimensioned zero is desired, supply it as the value for a constant polynomial; but if your vector type supports slicing, you'd better pass it as Polynomial({0: zero}) instead. The polynomial 1 + 2.x + 3.x**2 <- x may thus be generated by any of: Polynomial(1,2,3), Polynomial((1,2,3)), Polynomial([1,2,3]) or Polynomial({0:1,1:2,2:3}). The last, dictionary, form is generally preferrable for sparse polynomials (ie, rank >> the number of non-zero coefficients): tuple form is better for short and sweet forms like that illustrated. Note that setting z = Polynomial(0, 1) provides a `free variable' that can then be used to generate polynomials the way many folk prefer; with this, Polynomial(1,2,3) can simply be written 3*z*z +2*z +1. If you want a polynomial's representation to use the name z (rather than x), simply set the polynomial's .variablename to 'z' (but note that this only applies to that polynomial, not to ones computed from it).\n""" self.__coefs = {} # dictionary with natural keys try: arg, = args # do we have just one arg ? arg[:] # if so, is it a sequence ? self.__fromseq(arg) except ValueError: self.__fromseq(args) # several args except (AttributeError, TypeError, KeyError): # one non-sequence arg try: arg.items, arg.get(0, None) # is it a mapping ? except AttributeError: self.__store(0, arg) # constant polynomial else: self.__frombok(arg) # dictionary # self.__coefs should now be construed as immutable. # Coefficients only require the ability to add, multiply and divmod. def _get_coeff(val, oktypes=(types.IntType, types.ComplexType, types.FloatType, types.LongType)): if type(val) in oktypes: if val == long(val): return long(val) return val elif hasattr(val, '__add__') and hasattr(val, '__mul__'): # This means we can use polynomials as coefficients ... # likewise, linear maps, &c. return val else: raise invalidCoefficient def __store(self, power, coeff, g=_get_coeff): if coeff: self.__coefs[power] = g(coeff) else: self._zero = g(coeff) del _get_coeff def __fromseq(self, seq): i = 0 for v in seq: self.__store(i, v) i = 1 + i def __frombok(self, bok, ok=lambda k, i=(types.IntType, types.LongType): type(k) in i): for key, val in bok.items(): if not ok(key) or key < 0: raise unNaturalPower else: self.__store(key, val) def _lazy_get_rank_(self, ignored): for key in self.__coefs.keys(): if self.__coefs[key] == 0.0: del self.__coefs[key] try: return max(self.__coefs.keys()) except ValueError: return -1 def hcf(self, *others): for other in others: # Euclid's algorithm while other: self, other = other, self % other return self def coprime(self, *others): # coprime <=> hcf is a non-zero scalar return self.hcf(others).rank == 0 def _lazy_get_normalised_(self, ignored): if self.rank < 0: raise ValueError, "Can't normalise zero" scale = self.__coefs[self.rank] if scale == 1: return self return self / scale def __repr__(self): try: return self.__repr except AttributeError: pass names = [ self.variablename ] text = self.__represent(names, 0) # may grow names ! lamb = 'lambda %s: ' % ', '.join(names) if not text: ans = lamb + '0' else: ans = lamb + text self.__repr = ans return ans __str__ = __repr__ variablename = 'z' def format(num, names, depth): try: if num.imag == 0: num = num.real except AttributeError: pass try: return num.__represent(names, 1+depth) except AttributeError: return str(num) def __represent(self, names, depth, fmt=format): while depth >= len(names): if names[-1][0] == 'a': names.append('Z') else: names.append(chr(ord(names[-1][0]) - 1)) result, name = '', names[depth] for key in self._powers: val = self.coefficient(key) if key: if val == 1: frag = '' elif val == -1: frag = '-' else: frag = fmt(val, names, depth) if ' ' in frag: frag = '(' + frag + ')*' else: frag += '*' if key == 1: frag += name else: frag += '%s**%d' % (name, key) else: frag = fmt(val, names, depth) if frag[:1] == '-': result += ' ' + frag else: result += ' +' + frag if result[:2] == ' +': return result[2:] elif result[:1] == ' ': return result[1:] return result del format def __eachattr(self, each): bok = {} for k, v in self.__coefs.items(): bok[k] = each(v) return Polynomial(bok) def toreal(val): try: return val.real except AttributeError: return val def _lazy_get_real_(self, ignored, as=toreal): return self.__eachattr(as) del toreal def toimag(val): try: return val.imag except AttributeError: return 0 def _lazy_get_imag_(self, ignored, as=toimag): return self.__eachattr(as) del toimag def conjug8(val): try: return val.conjugate except AttributeError: return val def _lazy_get_conjugate_(self, ignored, as=conjug8): return self.__eachattr(as) del conjug8 def _lazy_get__powers_(self, ig): keys = self.__coefs.keys() keys.sort() keys.reverse() return tuple(keys) def __add__(self, other): try: sum = other.__coefs.copy() except AttributeError: sum = {0: other} zero = self._zero for k, v in self.__coefs.items(): sum[k] = sum.get(k, zero) + v return Polynomial(sum) __radd__ = __add__ def __sub__(self, other): sum, zero = self.__coefs.copy(), self._zero try: bok = other.__coefs except AttributeError: sum[0] = sum.get(0, zero) - other else: for key, val in bok.items(): sum[key] = sum.get(key, zero) - val return Polynomial(sum) def __rsub__(self, other): try: sum = other.__coefs.copy() except AttributeError: sum = {0: other} zero = self._zero for k, v in self.__coefs.items(): sum[k] = sum.get(k, zero) - v return Polynomial(sum) def __mul__(self, other): term = {} try: bok = other.__coefs except AttributeError: for key, val in self.__coefs.items(): term[key] = val * other else: zero = self._zero * other._zero for key, val in self.__coefs.items(): for cle, lue in bok.items(): sum = key + cle term[sum] = term.get(sum, zero) + val * lue return Polynomial(term) __rmul__ = __mul__ # abelian multiplication def __div__(self, other): q, r = self.__divmod__(other) if r and not r.__istiny(self._bigcoef): raise ValueError(self, 'not a multiple of', other, r) return q __floordiv__ = __div__ def __mod__(self, other): return self.__divmod__(other)[1] def divide(num, div): rat = num / div if div * rat == num: return rat return num * 1. / div def __divmod__(self, other, ratio=divide): """Solves self = q.other + r for r of rank < other.rank: returns (q, r) This depends on our coefficients forming a field. """ try: top = other.rank except AttributeError: top, other = 0, Polynomial({0: other}) if top < 0: raise ZeroDivisionError o = other.coefficient(top) if top == 0: bok = {} for k, v in self.__coefs.items(): bok[k] = v / o return Polynomial(bok), Polynomial(0) q, r = Polynomial(0.0), self got = self.rank # We now reduce the rank of r (by at least 1) at each iteration, by # shifting other*x**(got-top) times a scalar from r to q*other; thus, as # rank r is finite, it must eventually descend to 0. while got >= top: m = r.coefficient(got) while m: # take several slices off, in case of arithmetic error ... scale = Polynomial({ got - top: ratio(m, o) }) q, r = q + scale, r - scale * other n = r.coefficient(got) if abs(n) * 10 > abs(m): print "Wiping %g x**%d in Polynomial.__divmod__" % ( r.__coefs[got], got ) del r.__coefs[got] # => forcefully set to zero break m = n assert r.rank < got got = r.rank # assert r == self - q * other # tends to fail on small errors return q, r def scalarroot(val, exp, odd): assert val try: if val > 0: pass elif odd: return - ((-val)**exp) else: val = val + 0j # even root of -ve; coerce to complex except TypeError: pass # complex return val ** exp def __root(self, num, ratio=divide, root=scalarroot): # solve self = ans ** num assert num > 0 == self.rank % num top = self.rank / num ans = Polynomial({ top: root(self.coefficient(self.rank), 1./num, num%2)}) res = self - ans ** num while top > 0 and res: next = res.rank - (num-1)*ans.rank if next < top: top = next else: raise ValueError ans = ans + Polynomial({ top: ratio(res.coefficient(res.rank), num * ans.coefficient(ans.rank)**(num-1))}) res = self - ans ** num if res: raise ValueError return ans del divide, scalarroot def whole(val): try: return int(val) except OverflowError: return long(val) def coefficient(self, key, int=whole): val = self.__coefs.get(key, self._zero) try: i = int(val) if val == i: return i except (TypeError, AttributeError): pass return val def __pow__(self, other, int=whole, ok=lambda i, t=(types.IntType, types.LongType): type(i) in t): # Require other to be natural result, wer = 1, self try: if other.rank < 1: other = other.coefficient(0) # For some bizarre reason, x**0 isn't evaluated as x.__pow__(0) ! # When evaluating x ** 0 I find other == Polynomial(0) instead. except AttributeError: pass try: if other < 0: raise TypeError if self.rank < 0: return self # zero**+ve is zero m = 0 if not ok(other): i = int(other) if i == other: other = i else: # bad, unless we're very lucky i = other * self.rank if i != int(i): raise TypeError m = 1L except (AttributeError, TypeError, ValueError): raise unNaturalPower, other if m: # OK, we *might* be able to get away with this ... while 1: m = 1 + m i = m * other if i == int(i): break assert m <= self.rank # ... other is i/m; see if we have an m-th root: other, wer = int(i), self.__root(m) # will ValueError if not possible while other >= 1: if other % 2: result = result * wer wer, other = wer * wer, other / 2 return result del whole # use << as repeated integration, >> as repeated differentiation def __lshift__(self, other): """f << n -> nth integral of f""" 1L << other # evaluate to raise suitable error if any while other > 0: self, other = self.integral(), other - 1 return self def __rshift__(self, other): """f >> n -> nth derivative of f""" 1L >> other # evaluate in order to raise suitable error if any while other > 0: self, other = self.derivative, other - 1 return self def _lazy_get_derivative_(self, ignored): """Differentiate a polynomial. """ bok = {} for k, v in self.__coefs.items(): if k: bok[k-1] = k * v return Polynomial(bok) def integral(self, start=0, base=0): """Integrate a polynomial. Optional arguments, start and base, specify the constant of integration; the integral's value at start (whose default is 0) will be base (whose default is also zero). Note that self.integral(base=h) is equivalent to self.integral()+h; and self.integral(start, 0)(end) is the integral of self from start to end.""" bok = {} for k, v in self.__coefs.items(): bok[k+1] = v / (1.+k) ans = Polynomial(bok) # assert: ans(0) == 0 if start: ans = ans - ans(start) if base: ans = ans + base return ans # assert: self.integral().derivative == self # assert: self.derivative.integral(x, self(x)) == self, for any x # non-zero: safe and unambiguous def __nonzero__(self): return self.rank >= 0 # Comparison: of debatable value def __cmp__(self, other): return (self - other).sign or 0 # *Stronger* test for equality ... def __eq__(self, other): return (self - other).rank < 0 def __pos__(self): return self def __neg__(self): return 0 - self def __coerce__(self, other, boktyp=types.DictionaryType): try: if type(other.__coefs) == boktyp: return self, other except AttributeError: pass try: return self, Polynomial(other) except (unNaturalPower, invalidCoefficient): return None # Only useful if we want polys as keys in dictionaries ... def _lazy_get__lazy_hash_(self, ignored): result = 0 for key, val in self.__coefs.items(): result = result ^ hash(key) ^ hash(val) return result # How to evaluate a polynomial ... # For any T supporting # ({(T: :T)}: * :T), ({(V::T)}: * :{coefficients}) and ({(V: :V)}: + :V) # evaluation is ({(V: :T)}: :{polynomials}) # but the following should also suffice ... def __call__(self, arg): """Evaluate a polynomial as a function For a polynomial, p, with coefficients in some domain V, and a value t in some domain T supporting *: T-> V-> V and +: V-> V-> V, we can evaluate p(t) by substituting t in as the value of p's formal parameter.\n""" keys = self._powers try: top, keys = keys[0], keys[1:] except IndexError: return self._zero result = self.coefficient(top) for key in keys: # highest first while top > key: result, top = result * arg, top - 1 result = result + self.coefficient(key) while top > 0: result, top = result * arg, top - 1 return result # Calling one polynomial with another as input yields the composite of the two; # the following explores undoing that: def unafter(self, other): """Returns a polynomial which, if fed other, will yield self. The nature of polynomial arithmetic is such that, if f and g are polynomials, their composite (: g(f(x)) ←x :) is simply g(f). This method seeks to express self as g(other) for some g. Requires rank of self to be a multiple of rank of other, among other things; raises ValueError if the goal can't be met.\n""" residue, result = self, {} while residue: if residue.rank % other.rank: raise ValueError('Cannot factorise', residue, 'via', other) p = residue.rank / other.rank q, residue = divmod(residue, other ** p) assert q.rank == 0 assert residue.rank < p * other.rank result[p] = q.coefficient(0) return Polynomial(result) # unbefore would seem equivalent to arbitrary root-finding ! def _lazy_get_Gamma_(self, ignored): """Integrates self * (: exp(-t) ←t :{positives}). When self is power(n) the result is Gamma(n+1) = n!, so the result is just the sum of self's coefficients, each multiplied by the factorial of its order. The name is slighly misguided but not unreasonable. See atomic.py for a sub-class using this.\n""" n = self.rank ans = self.coefficient(n) while n > 0: ans, n = ans * n, n - 1 ans = ans + self.coefficient(n) return ans def _lazy_get_assquares_(self, ignored): """Decompose self as a sum of scaled squares, as far as possible. Returns a twople, (bok, poly) in which: poly is zero or a polynomial of odd rank; bok is a mapping from polynomials to scalars; if each key of bok is squared and multiplied by the corresponding value, summing the results and adding poly will yield self. """ bok, rem = {}, self while rem.rank % 2 == 0: i = rem.rank k, i = rem.coefficient(i), i / 2 assert k != 0 z = Polynomial({i: 1}) while i > 0: i = i - 1 q, ign = divmod(rem - k*z*z, 2*z*k) assert q.rank <= i z = z + Polynomial({i: q.coefficient(i)}) bok[z], rem = k, rem - k * z * z #print bok, rem assert self == rem + reduce(lambda y, (x, k): y + x*x*k, bok.items(), 0) return bok, rem def __pure_real(self): for v in self.__coefs.values(): try: if v.imag: return None except AttributeError: pass return 1 def _lazy_get__zero_(self, ig): try: return self.__coefs.values()[0] * 0 except IndexError: return 0 def _lazy_get_isreal_(self, ignored): return self.__pure_real() def _lazy_get_sign_(self, ignored): if self.rank < 0: return 0 # definitely everywhere zero if self.rank % 1 or not self.isreal: return None if self.rank == 0: return cmp(self.coefficient(0), 0) b, r = self.assquares if r.rank > 0: return None row = b.values() if r.rank == 0: row.append(r.coefficient(0)) lo, hi = min(row), max(row) if lo > 0: return +1 # definitely everywhere positive if hi < 0: return -1 # definitely everywhere negative return None # this is over-cautious from cardan import cubic def __cubic_root(self, cub=cubic): # for cubics, we know how to be exact ... ans = apply(cub, map(self.coefficient, (3, 2, 1, 0))) # Prefer real answers: nice = filter(lambda x: (x + 0j).imag == 0, ans) if nice: ans = nice # Prefer whole answers: nice = filter(lambda x: x == long(x), ans) if nice: ans = nice # Prefer bigger answers: hi, ans = ans[0], ans[1:] for lo in ans: if abs(lo) > abs(hi): hi = lo return hi del cubic def seek_root(self, guess=0, tol=1e-6): """Find an input which self maps to zero. Optional arguments: guess -- where to start searching [default: 0] tol -- tolerance [default: 1e-6]; if abs(self(x)) < tol, x is a root These are ignored if self is a quadratic or a real cubic; cardan.py's tools then provide for exact solution. Otherwise, Newton-Raphson is deployed. If self is everywhere or nowhere zero, ValueError is raised; this can only happen for constant polynomials. """ if self.rank < 1: raise ValueError, 'constant function is either nowhere or everywhere zero' if self.rank < 3 or (self.rank == 3 and self.__pure_real()): # Exploit exact solution: return self.__cubic_root() if self.coefficient(0) == 0: return 0 # easy special case to spot ;^) if self.normalised.sign and (guess + 0j).imag == 0: guess = guess + 1j # a pure real start point won't work grade, off = self.derivative, self(guess) while abs(off) > tol: d = grade(guess) #print 'Trying:', guess, '->', off, '@', d if d: guess = guess - off * 1. / d else: # at a stationary point try: guess = wobble + guess except NameError: if not self.isreal or self.sign: spin = 1 +1j else: spin = -2 #print 'Spinning:', spin wobble = spin guess = guess + wobble wobble = wobble * spin #print 'Wobbling:', wobble off = self(guess) return guess def seek_factor(self, guess=None): # do newton-raphson in the vector space of quadratic factors # stay within real if all coefficients are real if self.rank < 2: return self.normalised try: if guess.rank < 1: guess = None except AttributeError: guess = None if guess is None: try: root = self.seek_root() except ValueError: assert self.rank == 2 and self.isreal return self.normalised if self.isreal: try: if root.imag: guess = Polynomial({2: 1, 1: -2 * root.real, 0: root.real**2 +root.imag**2}) except AttributeError: pass if guess is None: guess = Polynomial({1: 1, 0: -root}) tiny = 1e-6 * self._bigcoef index = range(guess.rank) # not including rank itself; # leave guess' 1st coefficient alone; and r's rank is less than that of guess. while 1: r = self % guess if r.rank < 0: return guess # exact factor if r._bigcoef < tiny: return guess # pretty close small = 1e-3 * min(r._bigcoef, guess._bigcoef) matrix = [] for i in index: q = guess + Polynomial({i: small}) d = (r - self % q) / small row = [] for j in index: row.append(d.coefficient(j)) matrix.append(row) raise NotImplementedError, 'Need to divide vector by matrix' guess = guess - r / matrix # implement this ... def _lazy_get__bigcoef_(self, ignored): return max(map(abs, self.__coefs.values())) def __istiny(self, scale=1, maxrank=0): if self.rank > maxrank: return None if self.rank < 0 or self._bigcoef < scale * 1e-6: return 1 return None def _lazy_get_factors_(self, ignored): ans = [] while self.rank > 0: #print ans, self bok, rem = self.assquares if rem: if bok: f = apply(rem.hcf, bok.keys()) if f.rank > 0: f = f.normalised ans, self = ans + list(f.factors), self / f continue else: row = bok.keys() # assert: non-empty f = apply(row[0].hcf, row[1:]) if f.rank > 0: f = f.normalised ans, self = ans + list(2 * f.factors), self / f / f continue if len(bok) == 2: (X, a), (Y, b) = bok.items() if b < 0 < a: X, a, Y, b = Y, b, X, a if a < 0 < b: # difference of two squares b, a = b**.5, (-a)**.5 # self = (b.Y)**2 -(a.X)**2 X, Y = (b*Y -a*X).normalised, (b*Y +a*X).normalised ans, self = ans + list(X.factors + Y.factors), self / X / Y continue f = self.seek_factor() while 1: new, rem = divmod(self, f) if not rem.__istiny(self._bigcoef): break ans.append(f) self = new if self != 1: ans.append(self.coefficient(0)) return tuple(ans) del types, Lazy