sympycore/polynomials/algebra.py

   1 #
   2 # Author: Pearu Peterson
   3 # Created: February 2008
   4 #
   5 """ Provides PolynomialRing class.
   6 """
   7
   8 __docformat__ = "restructuredtext"
   9 __all__ = "PolynomialRing"
  10
  11 from ..core import classes, IntegerList, Expr
  12 from ..utils import SYMBOL, NUMBER, ADD, MUL, POW, SPARSE_POLY, SUB, DIV
  13 from ..basealgebra.algebra import Algebra
  14 from ..ring import CommutativeRing
  15 from ..basealgebra.verbatim import Verbatim
  16 from ..arithmetic.numbers import div
  17 from ..arithmetic.number_theory import multinomial_coefficients
  18
  19 def cmp_symbols(x, y):
  20     return cmp(str(x), str(y))
  21
  22 class PolynomialRingFactory(type):
  23     """ Factory of polynomial rings with symbols and coefficient ring.
  24     """
  25     def __new__(typ, name, bases, attrdict):
  26         if not attrdict.has_key('ring'):
  27             attrdict['ring'] = classes.Calculus
  28         if not attrdict.has_key('variables'):
  29             attrdict['variables'] = ()
  30             attrdict['nvars'] = 0
  31         cls = type.__new__(typ, name, bases, attrdict)
  32         cls.zero = cls.Number(0)
  33         cls.one = cls.Number(1)
  34         return cls
  35
  36     def __eq__(self, other):
  37         if isinstance(other, PolynomialRingFactory):
  38             return self.ring==other.ring and self.variables==other.variables
  39         return False
  40
  41     def __ne__(self, other):
  42         return not self==other
  43
  44     def __getitem__(self, ring_info, cache={}):
  45         """ Return a new polynomial ring class
  46
  47         Examples::
  48
  49           PolynomialRing[<seq of variables>, <coefficient ring>]
  50           PolynomialRing[<n>, <coefficient ring>] is
  51                 PolynomialRing[['X%s'%i for i in range(<n>)], <coefficient ring>]
  52           PolynomialRing[<variable info>] is
  53                 PolynomialRing[<variable info>, Calculus]
  54         """
  55         if isinstance(ring_info, (int, long)):
  56             nvars = ring_info
  57             assert nvars>=0,`nvars`
  58             variables = ['X%s'%i for i in range(nvars)]
  59             ring = classes.Calculus
  60         elif isinstance(ring_info, tuple) and isinstance(ring_info[-1], type):
  61             if len(ring_info)==2:
  62                 var_info, ring = ring_info
  63                 if isinstance(var_info,(int, long)):
  64                     nvars = var_info
  65                     assert nvars>=0,`nvars`
  66                     variables = ['X%s'%i for i in range(nvars)]
  67                 elif isinstance(var_info,(tuple, list, set)):
  68                     variables = var_info
  69                 elif isinstance(var_info, str):
  70                     variables = var_info,
  71                 else:
  72                     raise TypeError(`ring_info`)
  73             else:
  74                 variables = ring_info[:-1]
  75                 ring = ring_info[-1]
  76         elif isinstance(ring_info, (tuple, list, set)):
  77             variables = ring_info
  78             ring = classes.Calculus
  79         elif isinstance(ring_info, str):
  80             variables = ring_info,
  81             ring = classes.Calculus
  82         else:
  83             raise TypeError(`ring_info`)
  84         variables = tuple(sorted(variables, cmp=cmp_symbols))
  85         nvars = len(variables)
  86
  87         r = None #cache.get((variables, ring))
  88         if r is None:
  89             name = '%s[%s, %s]' % (self.__name__, tuple(variables), ring.__name__)
  90             r = PolynomialRingFactory(name,
  91                                       (self,),
  92                                       dict(nvars=nvars, ring = ring,
  93                                            variables = variables))
  94             #cache[variables, ring] = r
  95         return r
  96
  97     def is_subring(self, other):
  98         """ Check if self contains other as a subring, i.e. whether
  99         the other instances can be converted to self instances.
 100         """
 101         if not isinstance(other, PolynomialRingFactory):
 102             return False
 103         if self.ring != other.ring:
 104             return False
 105         return set(other.variables).issubset(self.variables)
 106
 107 AdditiveTuple = IntegerList
 108
 109 class PolynomialRing(CommutativeRing):
 110     """ Base class to polynomial rings that holds polynomial information
 111     using pairs ``(<exponents>: <coefficient>)`` stored in Python dictionary.
 112
 113     Suitable for representing sparse multivariate polynomials.
 114     """
 115
 116     __slots__ = ['_degree', '_ldegree']
 117     _degree = None
 118     _ldegree = None
 119     _str_value = None
 120
 121     __metaclass__ = PolynomialRingFactory
 122
 123     __str__ = CommutativeRing.__str__
 124     __repr__ = CommutativeRing.__repr__
 125
 126     @classmethod
 127     def convert(cls, data, typeerror=True):
 128         if isinstance(data, list):
 129             if cls.nvars==1 and data:
 130                 if type(data[0]) not in [tuple,list,set]:
 131                     data = [(i,c) for (i,c) in enumerate(data) if c]
 132             data = dict(data)
 133             return cls(SPARSE_POLY, data)
 134         elif isinstance(data, dict):
 135             return cls(SPARSE_POLY, data)
 136         r = super(CommutativeRing, cls).convert(data, typeerror=True)
 137         return r
 138
 139     @classmethod
 140     def get_operand_algebra(cls, head, index=0):
 141         if head in [ADD, SUB, MUL]:
 142             return cls
 143         if head is POW:
 144             if index==0:
 145                 return cls
 146             if index==1:
 147                 return int
 148         return cls.handle_get_operand_algebra_failure(head, index)
 149
 150     def __eq__(self, other):
 151         if type(other)==type(self):
 152             return self.pair == other.pair
 153         return CommutativeRing.__eq__(self, other)
 154
 155     @classmethod
 156     def Symbol(cls, obj):
 157         """ Return symbol element of a polynomial ring. The result
 158         may be an instance of a super polynomial ring.
 159
 160         Examples::
 161
 162           r = PolynomialRing['x']
 163           r.Symbol('x') -> r({1:1})
 164           r.Symbol('y') -> PolynomialRing['x','y']({(0,1):1})
 165         """
 166         try:
 167             i = list(cls.variables).index(obj)
 168         except ValueError:
 169             i = None
 170         if i is None:
 171             cls = PolynomialRing[cls.variables+(obj,), cls.ring]
 172             i = list(cls.variables).index(obj)
 173         l = [0]*cls.nvars
 174         l[i] = 1
 175         return cls(SPARSE_POLY, {AdditiveTuple(l):1})
 176
 177     @classmethod
 178     def convert_symbol(cls, obj):
 179         try:
 180             i = list(cls.variables).index(obj)
 181         except ValueError:
 182             i = None
 183         if i is None:
 184             cls = PolynomialRing[cls.variables+(obj,), cls.ring]
 185             i = list(cls.variables).index(obj)
 186         l = [0]*cls.nvars
 187         l[i] = 1
 188         return cls(SPARSE_POLY, {AdditiveTuple(l):1})
 189
 190     @classmethod
 191     def Number(cls, obj):
 192         """ Return number element of a polynomial ring.
 193
 194         Examples::
 195
 196           r = PolynomialRing['x']
 197           r.Number(2) -> r({0:2})
 198         """
 199         data = {AdditiveTuple((0,)*cls.nvars): obj} if obj else {}
 200         return cls(SPARSE_POLY, data)
 201
 202     @classmethod
 203     def Add(cls, *seq):
 204         r = cls(SPARSE_POLY, {})
 205         for t in seq:
 206             tcls = t.__class__
 207             if cls==tcls:
 208                 iadd_POLY_POLY(r, t, cls)
 209             elif cls.is_subring(tcls):
 210                 t = t.as_algebra(cls)
 211                 assert cls==t.__class__,`cls,t`
 212                 iadd_POLY_POLY(r, t, cls)
 213             elif tcls.is_subring(cls):
 214                 cls = tcls
 215                 r = r.as_algebra(cls)
 216                 iadd_POLY_POLY(r, t, cls)
 217             elif cls.ring==tcls.ring:
 218                 cls = PolynomialRing[cls.variables+tcls.variables, cls.ring]
 219                 r = r.as_algebra(cls)
 220                 t = t.as_algebra(cls)
 221                 iadd_POLY_POLY(r, t, cls)
 222             else:
 223                 raise NotImplementedError(`r, t`)
 224         if not r.data:
 225             return cls.Number(0)
 226         return r
 227
 228     @classmethod
 229     def Mul(cls, *seq):
 230         r = cls.one
 231         for t in seq:
 232             tcls = t.__class__
 233             if cls==tcls:
 234                 r = mul_POLY_POLY(r, t, cls)
 235             elif cls.is_subring(tcls):
 236                 t = t.as_algebra(cls)
 237                 assert cls==t.__class__,`cls,t`
 238                 r = mul_POLY_POLY(r, t, cls)
 239             elif tcls.is_subring(cls):
 240                 cls = tcls
 241                 r = r.as_algebra(cls)
 242                 r = mul_POLY_POLY(r, t, cls)
 243             elif cls.ring==tcls.ring:
 244                 cls = PolynomialRing[cls.variables+tcls.variables, cls.ring]
 245                 r = r.as_algebra(cls)
 246                 t = t.as_algebra(cls)
 247                 r = mul_POLY_POLY(r, t, cls)
 248             else:
 249                 raise NotImplementedError(`r, t`)
 250         if not r.data:
 251             return cls.Number(1)
 252         return r
 253
 254     @classmethod
 255     def Pow(cls, base, exp):
 256         r = cls.__pow__(base, exp)
 257         if r is NotImplemented:
 258             raise NotImplementedError(`base,exp`)
 259         return r
 260
 261     @classmethod
 262     def convert_coefficient(cls, obj, typeerror=True):
 263         if not isinstance(obj, cls.ring):
 264             r = cls.ring.convert_coefficient(obj, typeerror=False)
 265             if r is not None:
 266                 return r
 267             r = cls.ring.convert(obj, typeerror=False)
 268             if r is not None:
 269                 return r
 270         if typeerror:
 271             raise TypeError('%s.convert_coefficient: failed to convert %s instance'\
 272                             ' to coefficient algebra, expected int|long object'\
 273                             % (cls.__name__, obj.__class__.__name__))
 274         else:
 275             return NotImplemented
 276
 277     def as_verbatim(self):
 278         data = self.data
 279         symbols = [Verbatim.convert(v) for v in self.variables]
 280         terms = []
 281         for exps in sorted(data.keys(), reverse=True):
 282             coeff = data[exps]
 283             if coeff==1:
 284                 factors = []
 285             else:
 286                 factors = [Verbatim.convert(coeff)]
 287             if not isinstance(exps, tuple):
 288                 exps = exps,
 289             for s,e in zip(symbols, exps):
 290                 if e:
 291                     if e==1:
 292                         factors.append(s)
 293                     else:
 294                         factors.append(s**e) #XXX: ????
 295             if not factors:
 296                 terms.append(Verbatim.convert(coeff))
 297             elif len(factors)==1:
 298                 terms.append(factors[0])
 299             else:
 300                 terms.append(Verbatim(MUL,tuple(factors)))
 301         if not terms:
 302             return Verbatim.convert(0)
 303         if len(terms)==1:
 304             return terms[0]
 305         return Verbatim(ADD, tuple(terms))
 306
 307     def as_algebra(self, target_cls):
 308         cls = self.__class__
 309         if cls is target_cls:
 310             return self
 311         if issubclass(target_cls, PolynomialRing):
 312             target_ring = target_cls.ring
 313             if cls.variables == target_cls.variables:
 314                 data = dict([(e,target_ring.convert(c)) for e,c in self.data.iteritems()])
 315                 return target_cls(data)
 316             if cls.ring == target_ring:
 317                 new_cls = PolynomialRing[set(cls.variables + target_cls.variables), target_ring]
 318                 new_variables = list(new_cls.variables)
 319                 indices = [new_variables.index(v) for v in cls.variables]
 320                 data = {}
 321                 n = new_cls.nvars
 322                 for e,c in self.data.iteritems():
 323                     new_e = [0] * n
 324                     if not isinstance(e, tuple):
 325                         e = e,
 326                     for i,j in enumerate(indices):
 327                         new_e[j] = e[i]
 328                     data[AdditiveTuple(new_e)] = c
 329                 return new_cls(data)
 330         return NotImplemented
 331
 332     def __neg__(self):
 333         head, data = self.pair
 334         return type(self)(head, dict([(e,-c) for e, c in data.iteritems()]))
 335
 336     def __add__(self, other):
 337         cls = self.__class__
 338         other_cls = other.__class__
 339         if cls == other_cls and self.head is other.head:
 340             return add_POLY_POLY(self, other, cls)
 341         other = self.convert(other)
 342         if other is NotImplemented:
 343             return other
 344         other_cls = other.__class__
 345         if cls == other_cls and self.head is other.head:
 346             return add_POLY_POLY(self, other, cls)
 347         elif self.nvars < other.nvars:
 348             return other + self
 349         return self.head.add(cls, self, other)
 350
 351     __radd__ = __add__
 352
 353     def __mul__(self, other):
 354         cls = self.__class__
 355         other_cls = other.__class__
 356         if cls == other_cls and self.head is other.head:
 357             return mul_POLY_POLY(self, other, cls)
 358         other2 = self.convert_coefficient(other, typeerror=False)
 359         if other2 is not NotImplemented:
 360             return mul_POLY_COEFF(self, other2, cls)
 361         other = self.convert(other, typeerror=False)
 362         if other is NotImplemented:
 363             return other
 364         other_cls = other.__class__
 365         if cls == other_cls and self.head is other.head:
 366             return mul_POLY_POLY(self, other, cls)
 367         elif self.nvars < other.nvars:
 368             return other * self
 369         return self.head.commutative_mul(cls, self, other)
 370
 371     __rmul__ = __mul__
 372
 373     def __divmod__(self, other):
 374         cls = self.__class__
 375         other_cls = other.__class__
 376         if cls == other_cls:
 377             if cls.nvars==1:
 378                 return divmod_POLY1_POLY1_SPARSE(self, other, cls)
 379         return NotImplemented
 380
 381     def __div__(self, other):
 382         cls = self.__class__
 383         r = self.convert_coefficient(other, typeerror=False)
 384         if r is not NotImplemented:
 385             return mul_POLY_COEFF(self, div(1, r), cls)
 386         return divmod(self, other)[0]
 387
 388     def __mod__(self, other):
 389         return divmod(self, other)[1]
 390
 391     def __pow__(self, other):
 392         if isinstance(other, type(self)):
 393             if other.head is NUMBER:
 394                 other = other.data
 395             if other.head is SPARSE_POLY:
 396                 data = other.data
 397                 if len(data)==1:
 398                     exps, coeff = data.items()[0]
 399                     if exps.data==[0]*len(exps):
 400                         other = coeff
 401         if isinstance(other,(int, long)):
 402             return pow_POLY_INT(self, other, self.__class__)
 403         print `self`, `other.pair`
 404         return NotImplemented
 405
 406     def __call__(self, *args):
 407         # XXX: Use Horner scheme.
 408         r = 0
 409         if self.nvars==1:
 410             for exps, coeff in self.data.iteritems():
 411                 r = coeff * args[0]**exps + r
 412         else:
 413             for exps, coeff in self.data.iteritems():
 414                 r = reduce(lambda x,y: x*y, [a**e for a,e in zip(args,exps)], coeff) + r
 415         return r
 416
 417     @property
 418     def degree(self):
 419         d = self._degree
 420         if d is None:
 421             data = self.data
 422             if not data:
 423                 self._degree = d = 0
 424             elif self.nvars==0:
 425                 self._degree = d = 0
 426             elif self.nvars==1:
 427                 self._degree = d = max(data.keys())
 428             else:
 429                 self._degree = d = map(max,zip(*data.keys()))
 430         return d
 431
 432     @property
 433     def ldegree(self):
 434         d = self._ldegree
 435         if d is None:
 436             data = self.data
 437             if not data:
 438                 self._ldegree = d = 0
 439             elif self.nvars==0:
 440                 self._ldegree = d = 0
 441             elif self.nvars==1:
 442                 self._ldegree = d = min(data.keys())
 443             else:
 444                 self._ldegree = d = map(min,zip(*data.keys()))
 445         return d
 446
 447     @property
 448     def coeff(self):
 449         data = self.data
 450         nvars = self.nvars
 451         if nvars==1:
 452             if data:
 453                 return data[self.degree]
 454             return self.zero
 455         raise NotImplementedError(`self,nvars`)
 456
 457     def diff(self, index=0): # TODO: don't use index as the order of variables is automatically sorted
 458         if self.nvars==0:
 459             return self.zero
 460         head, data = self.pair
 461         d = {}
 462         if self.nvars==1:
 463             assert index==0,`index`
 464             for exps, coeff in data.iteritems():
 465                 if exps:
 466                     d[exps-1] = coeff * exps
 467         else:
 468             assert 0<=index<self.nvars, `index`
 469             for exps, coeff in data.iteritems():
 470                 e = exps[index]
 471                 if e:
 472                     d[exps.add(index,-1)] = coeff * e
 473         return type(self)(head, d)
 474
 475     def variable_diff(self, variable, order=1):
 476         try:
 477             index = list(self.variables).index(variable)
 478         except ValueError:
 479             index = None
 480         if index is not None:
 481             return self.head.diff_index(type(self), self.data, self, index, order=order)
 482         raise NotImplementedError(`self.variables, variable, index`)
 483
 484     def variable_integrate(self, variable, *bounds):
 485         """ Return integral over variable.
 486
 487         poly.variable_integrate(x) -> indefinite integral
 488         poly.variable_integrate(x, a, b) -> definite integral
 489         """
 490         try:
 491             index = list(self.variables).index(variable)
 492         except ValueError:
 493             index = None
 494         if index is not None:
 495             indef_integral = self.head.integrate_indefinite_index(type(self), self.data, self, index)
 496             if bounds:
 497                 low, high = bounds
 498                 return indef_integral.variable_subs(variable, high) - indef_integral.variable_subs(variable, low)
 499             return indef_integral
 500         raise NotImplementedError(`self.variables, variable, index`)
 501
 502     def variable_subs(self, variable, newexpr):
 503         """ Substitute variable with newexpr that must be instance of the same ring.
 504         """
 505         cls = type(self)
 506         newexpr = cls(newexpr)
 507         try:
 508             index = list(self.variables).index(variable)
 509         except ValueError:
 510             index = None
 511         if index is not None:
 512             head, data = self.pair
 513             result = cls.Number(0)
 514             variables = cls.variables
 515             for exps, coeff in data.iteritems():
 516                 term = cls.Number(1)
 517                 for i,exp in enumerate(exps):
 518                     if exp:
 519                         if i==index:
 520                             term *= newexpr**exp
 521                         else:
 522                             term *= cls.Symbol(variables[i])**exp
 523                 result += term * cls.Number(coeff)
 524             return result
 525         raise NotImplementedError(`self.variables, variable, index`)
 526
 527 def divmod_POLY1_POLY1_SPARSE(lhs, rhs, cls):
 528     d2 = rhs.degree
 529     c2 = rhs.coeff
 530     if not c2:
 531         raise ZeroDivisionError, "polynomial division"
 532     other_iter = rhs.data.iteritems
 533     dq = {}
 534     dr = dict(lhs.data)
 535     dseq = [0] + dr.keys()
 536     dr_get = dr.get
 537     dseq_append = dseq.append
 538     dseq_remove = dseq.remove
 539     while 1:
 540         d1 = max(dseq)
 541         if d1 < d2:
 542             return cls(SPARSE_POLY, dq), cls(SPARSE_POLY, dr)
 543         c1 = dr[d1]
 544         e = d1 - d2
 545         c = div(c1, c2)
 546         dq[e] = c
 547         nc = -c
 548         for e0,c0 in other_iter():
 549             e0 = e0 + e
 550             c0 = nc * c0
 551             b = dr_get(e0)
 552             if b is None:
 553                 dr[e0] = c0
 554                 dseq_append(e0)
 555             else:
 556                 b = b + c0
 557                 if b:
 558                     dr[e0] = b
 559                 else:
 560                     del dr[e0]
 561                     dseq_remove(e0)
 562
 563 def divmod_POLY1_POLY1_DENSE(lhs, rhs, cls):
 564     if not rhs.coeff:
 565         raise ZeroDivisionError, "polynomial division"
 566     n = lhs.degree
 567     nv = rhs.degree
 568     u = [0]*(n+1)
 569     v = [0]*(nv+1)
 570     for e,c in lhs.data.iteritems():
 571         u[e] = c
 572     for e,c in rhs.data.iteritems():
 573         v[e] = c
 574     r = list(u)
 575     q = [0] * (len(r)+1)
 576     for k in range(n-nv, -1, -1):
 577         q[k] = div(r[nv+k], v[nv])
 578         for j in range(nv+k-1, k-1, -1):
 579             r[j] -= q[k]*v[j-k]
 580     for j in range(nv, n+1, 1):
 581         r[j] = 0
 582     q, r = cls(q), cls(r)
 583     return q, r
 584
 585 def pow_POLY_INT(base, exp, cls):
 586     if exp==0:
 587         return cls.one
 588     if exp==1:
 589         return base
 590     if exp < 0:
 591         return NotImplemented
 592     data = base.data
 593     if len(data)==1:
 594         exps, coeff = data.items()[0]
 595         return cls(SPARSE_POLY, {exps * exp: coeff ** exp})
 596     nvars = cls.nvars
 597     d = {}
 598     items = data.items()
 599     m = len(data)
 600     for k,c_kn in multinomial_coefficients(m, exp).iteritems():
 601         new_exps = AdditiveTuple((0,)*nvars)
 602         new_coeff = c_kn
 603         for i,e in enumerate(k):
 604             if e:
 605                 exps, coeff = items[i]
 606                 new_exps += exps * e
 607                 new_coeff *= coeff ** e
 608         b = d.get(new_exps)
 609         if b is None:
 610             d[new_exps] = new_coeff
 611         else:
 612             c = b + new_coeff
 613             if c:
 614                 d[new_exps] = c
 615             else:
 616                 del d[new_exps]
 617     return cls(SPARSE_POLY, d)
 618
 619
 620 def add_POLY_POLY(lhs, rhs, cls):
 621     d = dict(lhs.data)
 622     for exps, coeff in rhs.data.iteritems():
 623         b = d.get(exps)
 624         if b is None:
 625             d[exps] = coeff
 626         else:
 627             c = b + coeff
 628             if c:
 629                 d[exps] = c
 630             else:
 631                 del d[exps]
 632     return cls(SPARSE_POLY, d)
 633
 634 def iadd_POLY_POLY(lhs, rhs, cls):
 635     d = lhs.data
 636     for exps, coeff in rhs.data.iteritems():
 637         b = d.get(exps)
 638         if b is None:
 639             d[exps] = coeff
 640         else:
 641             c = b + coeff
 642             if c:
 643                 d[exps] = c
 644             else:
 645                 del d[exps]
 646
 647
 648 def mul_POLY_POLY(lhs, rhs, cls):
 649     d = {}
 650     for exps1, coeff1 in lhs.data.iteritems():
 651         for exps2, coeff2 in rhs.data.iteritems():
 652             exps = exps1 + exps2
 653             coeff = coeff1 * coeff2
 654             b = d.get(exps)
 655             if b is None:
 656                 d[exps] = coeff
 657             else:
 658                 c = b + coeff
 659                 if c:
 660                     d[exps] = c
 661                 else:
 662                     del d[exps]
 663     return cls(SPARSE_POLY, d)
 664
 665 def mul_POLY_COEFF(lhs, rhs, cls):
 666     d = {}
 667     for exps, coeff in lhs.data.iteritems():
 668         b = d.get(exps)
 669         if b is None:
 670             d[exps] = coeff * rhs
 671         else:
 672             c = b + coeff * rhs
 673             if c:
 674                 d[exps] = c
 675             else:
 676                 del d[exps]
 677     return cls(SPARSE_POLY, d)