research/directadd5.py

   1 import gmpy
   2
   3 NUMBER = 'N'
   4 SYMBOLIC = 'S'
   5 TERMS = '+'
   6 FACTORS = '*'
   7
   8 class rational(tuple):
   9
  10     def __new__(cls, p, q, tnew=tuple.__new__):
  11         x, y = p, q
  12         while y:
  13             x, y = y, x % y
  14         if x != 1:
  15             p //= x
  16             q //= x
  17         if q == 1:
  18             return p
  19         return tnew(cls, (p, q))
  20
  21     def __str__(self):
  22         return "%i/%i" % self
  23
  24     __repr__ = __str__
  25
  26     # not needed when __new__ normalizes to ints
  27     # __nonzero__
  28     # __eq__
  29     # __hash__
  30
  31     def __add__(self, other):
  32         p, q = self
  33         if isinstance(other, int):
  34             r, s = other, 1
  35         else:
  36             r, s = other
  37         return rational(p*s + q*r, q*s)
  38
  39     __radd__ = __add__
  40
  41     def __mul__(self, other):
  42         p, q = self
  43         if isinstance(other, int):
  44             r, s = other, 1
  45         elif isinstance(other, rational):
  46             r, s = other
  47         else:
  48             return NotImplemented
  49         return rational(p*r, q*s)
  50
  51     __rmul__ = __mul__
  52
  53
  54
  55 def make_int(n):
  56     return (NUMBER, n)
  57
  58 def make_rational(p, q):
  59     return (NUMBER, rational(p, q))
  60
  61 zero = make_int(0)
  62 one = make_int(1)
  63 two = make_int(2)
  64 half = make_rational(1, 2)
  65
  66 def make_symbol(name):
  67     return (SYMBOLIC, name)
  68
  69 def convert_terms(x):
  70     if x[0] is NUMBER:
  71         return (TERMS, [(one, x[1])])
  72     else:
  73         return (TERMS, [(x, 1)])
  74
  75 def add_default(a, b):
  76     if a[0] is not TERMS: a = convert_terms(a)
  77     if b[0] is not TERMS: b = convert_terms(b)
  78     return add_terms_terms(a, b)
  79
  80 def add_terms_terms(a, b):
  81     c = dict(a[1])
  82     for e, x in b[1]:
  83         if e in c:
  84             c[e] += x
  85             if not c[e]:
  86                 del c[e]
  87         else:
  88             c[e] = x
  89     if not c:
  90         return zero
  91     citems = c.items()
  92     if len(citems) == 1:
  93         w = citems[0]
  94         if w[1] == 1:
  95             return w[0]
  96     return (TERMS, frozenset(citems))
  97
  98 def add_symbolic_symbolic(a, b):
  99     if a == b:
 100         return (TERMS, frozenset([(a, 2)]))
 101     return (TERMS, frozenset([(a, 1), (b, 1)]))
 102
 103 def add_rational_rational(a, b):
 104     return (NUMBER, a[1]+b[1])
 105
 106 addition = {
 107     (TERMS, TERMS)       : add_terms_terms,
 108     (SYMBOLIC, SYMBOLIC) : add_symbolic_symbolic,
 109     (NUMBER, NUMBER) : add_rational_rational
 110 }
 111
 112
 113 def mul_default(a, b):
 114     if a[0] is not FACTORS: a = (FACTORS, [(a, 1)])
 115     if b[0] is not FACTORS: b = (FACTORS, [(b, 1)])
 116     return mul_factors_factors(a, b)
 117
 118 def mul_factors_factors(a, b):
 119     c = dict(a[1])
 120     for e, x in b[1]:
 121         if e in c:
 122             c[e] += x
 123             if not c[e]:
 124                 del c[e]
 125         else:
 126             c[e] = x
 127     if not c:
 128         return one
 129     return (FACTORS, frozenset(c.items()))
 130
 131 def mul_rational_symbolic(a, b):
 132     if a == zero: return zero
 133     if a == one: return b
 134     return (TERMS, frozenset([(b, a[1])]))
 135
 136 def mul_symbolic_rational(a, b):
 137     return mul_rational_symbolic(b, a)
 138
 139 def mul_rational_terms(a, b):
 140     if a == zero: return zero
 141     if a == one: return b
 142     p = a[1]
 143     return (TERMS, frozenset([(e, p*x) for e, x in b[1]]))
 144
 145 def mul_terms_rational(a, b):
 146     return mul_rational_terms(b, a)
 147
 148 def mul_rational_rational(a, b):
 149     return (NUMBER, a[1]*b[1])
 150
 151 multiplication = {
 152     (NUMBER, SYMBOLIC) : mul_rational_symbolic,
 153     (SYMBOLIC, NUMBER) : mul_symbolic_rational,
 154     (NUMBER, TERMS)    : mul_rational_terms,
 155     (TERMS, NUMBER)    : mul_terms_rational,
 156     (FACTORS, FACTORS)   : mul_factors_factors,
 157     (NUMBER, NUMBER) : mul_rational_rational
 158 }
 159
 160 def show(x):
 161     if isinstance(x, int):
 162         return str(x)
 163     h = x[0]
 164     if h is NUMBER:
 165         return str(x[1])
 166     if h is SYMBOLIC:
 167         return x[1]
 168     if h is TERMS:
 169         return ' + '.join((show(b)+"*"+show(a)) for (a, b) in x[1])
 170     if h is FACTORS:
 171         return ' * '.join(("(%s)**(%i/%i)" % (show(a), show(b))) for (a, b) in x[1])
 172     return str(x)
 173
 174 class Expr:
 175
 176     def __init__(self, value=None, _s=None):
 177         if _s:
 178             self._s = _s
 179             return
 180         if isinstance(value, int):
 181             self._s = (NUMBER, value)
 182         elif isinstance(value, str):
 183             self._s = make_symbol(value)
 184
 185     def __repr__(self):
 186         return show(self._s)
 187
 188     def __add__(self, other):
 189         if not isinstance(other, Expr):
 190             other = Expr(other)
 191         s, t = self._s, other._s
 192         add = addition.get((s[0], t[0]), add_default)
 193         return Expr(_s=add(self._s, other._s))
 194
 195     __radd__ = __add__
 196
 197     def __mul__(self, other):
 198         if not isinstance(other, Expr):
 199             other = Expr(other)
 200         s, t = self._s, other._s
 201         mul = multiplication.get((s[0], t[0]), mul_default)
 202         return Expr(_s=mul(s, t))
 203
 204     __rmul__ = __mul__
 205
 206 def Rational(p, q):
 207     return Expr(_s=make_rational(p, q))
 208
 209
 210 # use time() instead on unix
 211 import sys
 212 if sys.platform=='win32':
 213     from time import clock
 214 else:
 215     from time import time as clock
 216
 217 import sympycore as sympy
 218
 219 def time1():
 220     x = Expr('x')
 221     y = Expr('y')
 222     z = Expr('z')
 223     a = Rational(1,2)
 224     b = Rational(3,4)
 225     c = Rational(5,6)
 226     A = a*x + b*y + c*z
 227     B = b*x + c*y + a*z
 228     t1 = clock()
 229     n = 1000
 230     while n:
 231         3*(a*x+b*y+c*z)
 232         #A + B
 233         #x + y
 234         #a + b
 235         #a * b
 236         n -= 1
 237     t2 = clock()
 238     return 100 / (t2-t1)
 239
 240 def time2(n=1000):
 241     Symbol = sympy.Symbol
 242     Number = sympy.Number
 243     x = Symbol('x')
 244     y = Symbol('y')
 245     z = Symbol('z')
 246     a = Number(1,2)
 247     b = Number(3,4)
 248     c = Number(4,5)
 249     A = a*x + b*y + c*z
 250     B = b*x + c*y + a*z
 251     t1 = clock()
 252     while n:
 253         3*(a*x+b*y+c*z)
 254         #A + B
 255         #x + y
 256         #a + b
 257         #a * b
 258         n -= 1
 259     t2 = clock()
 260     return 100 / (t2-t1)
 261
 262
 263 def timing():
 264     t1 = time1()
 265     t2 = time2()
 266     return t1, t2, t1/t2
 267
 268 print "without psyco"
 269 print timing()
 270 print timing()
 271 print timing()
 272
 273 from sympycore import profile_expr
 274
 275 profile_expr('time2(1000)')
 276
 277 try:
 278     import psyco
 279 except ImportError:
 280     psyco = None
 281 if psyco is not None:
 282     psyco.full()
 283
 284     print "with psyco"
 285     print timing()
 286     print timing()
 287     print timing()