attachment:rewrite-20190430.sage of symbolics/rewrite

   1 """
   2 François Maltey - october 2010
   3 
   4 The "mapexpression" function maps recursively an expression as a tree.
   5 
   6 This function:
   7 
   8 - remains the framework of the expression,
   9 - treats first the leaves and then the composed expression, and
  10 - can change only some subtree.
  11 
  12 This change can be done:
  13 
  14 - at the main level only,
  15 - at some levels in the tree,
  16 - at all the levels.
  17 
  18 The function goes thru sum and product and can consider
  19 that the depth of the subexpression doesn't change.
  20 
  21 expr = the current expression
  22 mulDepth = 0 or 1, if 0 the depth in the tree remains the same in a sum.
  23 addDepth = 0 or 1, if 1 the depth in the tree increases in a product.
  24 fct = the effective function called for subtrees
  25 param = the parameter of fct
  26 level = -1 for a fully recursive call, or the list of levels to treat.
  27 level = [0] for the first main level.
  28 
  29 EXAMPLES::
  30 
  31     sage: rewrite(exp(x), "exp2sincos")
  32     -I*sin(I*x) + cos(I*x)
  33     sage: rewrite(exp(-I*x), "exp2sincos")
  34     -I*sin(x) + cos(x)
  35     sage: rewrite(exp(a+I*b), "exp2trig")
  36     (sinh(a) + cosh(a))*(I*sin(b) + cos(b))
  37     sage: rewrite((e^x)^2-e^(-2*x)+e^(-4*x)+(e^x)^4, 'exp2sinhcosh')
  38     2*sinh(2*x) + 2*cosh(4*x)
  39 
  40     sage: (binomial(n, k)*factorial(n)).rewrite(target=gamma)
  41     gamma(n+1)^2/(gamma(k+1)*gamma(n-k+1))
  42 
  43     sage: (exp(x)*tan(x)).rewrite(source=tan, target=sin)
  44     exp(x)*sin(x)*cos(x)
  45 """
  46 import operator
  47 
  48 from sage.functions.trig import (cos, sin, tan, cot, arctan, arcsin, arccos,
  49                                  sec, csc, arcsec, arccsc, arccot)
  50 from sage.functions.hyperbolic import (cosh, sinh, tanh, arctanh, arcsech,
  51                                        arccsch,  arccoth, arccosh, arcsinh,
  52                                        coth, csch, sech)
  53 from sage.functions.log import exp, log, ln
  54 from sage.functions.gamma import gamma
  55 from sage.functions.other import sqrt, factorial, floor
  56 from sage.rings.real_mpfr import RR
  57 from sage.symbolic.ring import SR
  58 from sage.misc.misc_c import prod
  59 from sage.all import pi, I
  60 
  61 
  62 def mapexpression(expr, fct, param, level, addDepth=0, mulDepth=0):
  63     def mapex(expr, depth):               # a very local function
  64         opor = expr.operator()
  65         opands = expr.operands()
  66         if (opor is None):
  67             return expr       # a leaf in the expression tree
  68         if (opor == operator.add):           # recursive call thru sum
  69             opands = map(lambda ex: mapex(ex, depth + addDepth), opands)
  70             return sum(opands)
  71         if (opor == operator.mul):           # recursive call thru mul
  72             opands = map(lambda ex: mapex(ex, depth + mulDepth), opands)
  73             return prod(opands)
  74         if (level == -1) or (level[-1] >= depth):  # recursive call over operands
  75             opands = map(lambda ex: mapex(ex, depth + 1), opands)
  76         if level == -1 or depth in level:  # root of the subtree must be changed
  77             return fct(opor, opands, param)
  78         return opor(*opands)  # opands may or may not be changed by a recursive call
  79     return mapex(expr, 0)
  80 
  81 
  82 def pseudoPositive(expr):
  83     """
  84     Choose the sign of expressions in respect of the opposite form.
  85 
  86     If a-b is said "positive" then -(a-b) is negative.
  87     """
  88     if expr._is_real():
  89         return bool(RR(expr) >= 0)  # can be improved by ._is_positive()
  90     if expr._is_numeric():      # see ._is_positive()
  91         return bool((expr.real() > 0) or (expr.real() == 0 and expr.imag() > 0))
  92     if expr._is_symbol():
  93         return True  # a_variable as x or a "is positive"
  94     opor = expr.operator()
  95     opands = expr.operands()
  96     if opor == operator.mul:            # look at the last factor in a product
  97         return pseudoPositive(opands[-1])  # it's the number if there is one
  98     if opor == operator.add:    # look at the first term in a sum
  99         return pseudoPositive(opands[0])
 100     return True    # for functions call as sin(x)...
 101 
 102 
 103 def pseudoRealAndImag(expr):
 104     """try to get (a,b) from a+i*b"""
 105     opands = expr.operands()     # don't decompose (a+i*b)*(x+i*y)
 106     opor = expr.operator()
 107     if opor == operator.mul:    # but treat Complex_number * expressions
 108         coeff = 1
 109         unit = 1
 110         for ex in opands:
 111             if ex._is_numeric():
 112                 coeff = coeff * ex
 113             else:
 114                 unit = unit * ex
 115         return (coeff.real() * unit, coeff.imag() * unit)
 116     elif opor == operator.add:  # and treat sum
 117         rp = 0
 118         ip = 0
 119         for ex in opands:
 120             if ex._is_numeric():
 121                 rp += ex.real()
 122                 ip += ex.imag()
 123             res = pseudoRealAndImag(ex)
 124             rp += res[0]
 125             ip += res[1]
 126         return (rp, ip)
 127     return (expr, 0)             # in doubt remain the expression as real
 128 
 129 
 130 def searchRule(source, target):
 131     """
 132     There are 2 uses of rewrite:
 133 
 134     - either by the 2 parameters source=... and target=...
 135     - either by the name "source2target" of this rule
 136 
 137     This function transforms the first method to the second one
 138     """
 139     if source in [sin, cos] and target in [sinh, cosh]:
 140         return ["trigo2trigh"]
 141     elif source == exp and target in [sin, cos]:
 142         return ["exp2trigo"]
 143     elif source == exp and target in [sinh, cosh]:
 144         return ["exp2trigh"]
 145     elif source in [sin, cos] and target == exp:
 146         return ["exp2trigo"]
 147     elif source in [sinh, cosh] and target == exp:
 148         return ["exp2trigh"]
 149     else:
 150         return "cannot match source and target, use one rule or a list of rules"
 151 
 152 
 153 def rewrite(expr, rule=None, source=None, target=None, filter=None, level=-1):
 154     if (rule is None) == (target is None):
 155         raise ValueError("must choose either rule=... or target=...")
 156     if target is None and source is not None:
 157         raise ValueError("must choose target=... for a defined source")
 158     if rule is None:
 159         rule = searchRule(source, target)
 160     elif type(rule) == str:
 161         rule = [rule]
 162     else:
 163         raise ValueError('a rule must be a string or a list of strings')
 164     if filter is None:
 165         def filter(ex):
 166             return True
 167     if type(filter).__name__ != 'function':
 168         return 'filter must be a function'
 169     if level != -1 and type(level) != list:
 170         raise ValueError('level must be -1 for a fully recursive '
 171                          'call or a list of integer')
 172     k = 0
 173     while k < len(rule):
 174         if rule[k] == "sinhcosh2exp":   # sinh(x) -> (exp(x)-exp(-x))/2, ...
 175             expr = mapexpression(expr, sinhcosh2exp, filter, level)
 176         elif rule[k] == "sincos2exp":   # sin(x) -> (exp(ix)-exp(-ix))/(2i), ...
 177             expr = mapexpression(expr, sincos2exp, filter, level)
 178         elif rule[k] == "exp2sinhcosh":  # exp(x) -> sinh(x)+cosh(x)
 179             expr = mapexpression(expr, exp2sinhcosh, filter, level)
 180         elif rule[k] == "exp2sincos":  # exp(x)=exp(i*(-ix)) -> cos(ix)-i*sin(ix)
 181             expr = mapexpression(expr, exp2sincos, filter, level)
 182 
 183         elif rule[k] == "trigo2trigh":  # cos(x) -> cosh(ix), ...
 184             expr = mapexpression(expr, trigo2trigh, filter, level)
 185         elif rule[k] == "trigh2trigo":  # cosh(x) -> cos(i*x), ...
 186             expr = mapexpression(expr, trigh2trigo, filter, level)
 187 
 188         elif rule[k] == "trigo2sincos":  # tan, cot, sec, csc -> sin, cos
 189             expr = mapexpression(expr, trigo2sincos, filter, level)
 190         elif rule[k] == "trigh2sinhcosh":  # tanh, coth, sech, csch -> sinh, cosh
 191             expr = mapexpression(expr, trigh2sinhcosh, filter, level)
 192 
 193         elif rule[k] == "exp2trig":
 194             # exp(a+i*b)->(cosh(a)+sinh(a))*(cos(b)+i*sin(b))
 195             expr = mapexpression(expr, exp2trig, filter, level)
 196         elif rule[k] == "lessIinExp":  # exp(a+i*b) -> exp(a)*(cos(b)+i*sin(b))
 197             expr = mapexpression(expr, lessIinExp, filter, level)
 198         elif rule[k] == "lessIinTrig":   # cos(i*x) -> cosh(x), ...
 199             expr = mapexpression(expr, lessIinTrig, filter, level)
 200 
 201         elif rule[k] == "trigo2exp":
 202             rule.extend(["sincos2exp", "trigo2sincos"])
 203         elif rule[k] == "trigh2exp":
 204             rule.extend(["sinhcosh2exp", "trigh2sinhcosh"])
 205         elif rule[k] == "trig2exp":
 206             rule.extend(["trigo2exp", "trigh2exp"])
 207 
 208         elif rule[k] == "cos22sin":
 209             expr = mapexpression(expr, squareInPow,
 210                                  [cos, filter, positiveCos,
 211                                   lambda ex: 1 - positiveSin(ex)**2], level)
 212         elif rule[k] == "sin22cos":
 213             expr = mapexpression(expr, squareInPow,
 214                                  [sin, filter, positiveSin,
 215                                   lambda ex:1-positiveCos(ex)**2], level)
 216         elif rule[k] == "cosh22sinh":
 217             expr = mapexpression(expr, squareInPow,
 218                                  [cosh, filter, positiveCosh,
 219                                   lambda ex: 1 + positiveSinh(ex)**2], level)
 220         elif rule[k] == "sinh22cosh":
 221             expr = mapexpression(expr, squareInPow,
 222                                  [sinh, filter, positiveSinh,
 223                                   lambda ex: positiveCosh(ex)**2-1], level)
 224         elif rule[k] == "tan22cos":
 225             expr = mapexpression(expr, squareInPow,
 226                                  [tan, filter, positiveTan,
 227                                   lambda ex: 1 / positiveCos(ex)**2-1], level)
 228         elif rule[k] == "cot22sin":
 229             expr = mapexpression(expr, squareInPow,
 230                                  [cot, filter, positiveCot,
 231                                   lambda ex: 1 / positiveSin(ex)**2-1], level)
 232         elif rule[k] == "tanh22cosh":
 233             expr = mapexpression(expr, squareInPow,
 234                                  [tanh, filter, positiveTanh,
 235                                   lambda ex: 1 - 1 / positiveCosh(ex)**2],
 236                                  level)
 237         elif rule[k] == "coth22sinh":
 238             expr = mapexpression(expr, squareInPow,
 239                                  [coth, filter, positiveCoth,
 240                                   lambda ex: 1 / positiveSinh(ex)**2-1], level)
 241         elif rule[k] == "cos22tan":
 242             expr = mapexpression(expr, squareInPow,
 243                                  [cos, filter, positiveCos,
 244                                   lambda ex: 1 / (positiveTan(ex)**2+1)], level)
 245         elif rule[k] == "tancot22sincos":
 246             rule.extend(["tan22cos", "cot22sin"])
 247         elif rule[k] == "tanhcoth22sinhcosh":
 248             rule.extend(["tanh22cosh", "coth22sinh"])
 249         elif rule[k] == "sin2tancos":
 250             expr = mapexpression(expr, sin2tancos, filter, level)
 251         elif rule[k] == "sincos2tan":
 252             expr = mapexpression(expr, sin2tancos, filter, level)
 253             expr = expr.rational_simplify()
 254             expr = mapexpression(expr, squareInPow,
 255                                  [cos, filter, positiveCos,
 256                                   lambda ex:1/(positiveTan(ex)**2+1)], level)
 257         elif rule[k] == "sinh2tanhcosh":
 258             expr = mapexpression(expr, sinh2tanhcosh, filter, level)
 259 
 260         elif rule[k] == "sincos2tanHalf":
 261             expr = mapexpression(expr, sincos2tanHalf, filter, level)
 262         elif rule[k] == "sinhcosh2tanhHalf":
 263             expr = mapexpression(expr, sinhcosh2tanhHalf, filter, level)
 264 
 265         elif rule[k] == "exp2pow":
 266             expr = mapexpression(expr, exp2pow, filter, level)
 267         elif rule[k] == "pow2exp":
 268             expr = mapexpression(expr, pow2exp, filter, level)
 269 
 270         elif rule[k] == "arcsin2arccos":
 271             expr = mapexpression(expr, arcsin2arccos, filter, level)
 272         elif rule[k] == "arccos2arcsin":
 273             expr = mapexpression(expr, arccos2arcsin, filter, level)
 274 
 275         elif rule[k] == "arctrigo2log":
 276             expr = mapexpression(expr, arctrigo2log, filter, level)
 277         elif rule[k] == "arctrigh2log":
 278             expr = mapexpression(expr, arctrigh2log, filter, level)
 279         elif rule[k] == "log2arctan":
 280             expr = mapexpression(expr, log2arctan, filter, level)
 281         elif rule[k] == "log2arctanh":
 282             expr = mapexpression(expr, log2arctanh, filter, level)
 283 
 284         elif rule[k] == "factorial2gamma":
 285             expr = mapexpression(expr, factorial2gamma, filter, level)
 286         elif rule[k] == "gamma2factorial":
 287             expr = mapexpression(expr, gamma2factorial, filter, level)
 288         elif rule[k] == "pseudoParity":
 289             expr = mapexpression(expr, pseudoParity, filter, level)
 290         elif rule[k] == "logMainBranch":
 291             expr = mapexpression(expr, logMainBranch, filter, level)
 292         elif rule[k] == "arctrigoMainBranch":
 293             expr = mapexpression(expr, arctrigoMainBranch, filter, level)
 294         else:
 295             return "unknown rule"
 296         k += 1
 297     return expr
 298 
 299 
 300 def positiveSin(expr):
 301     if pseudoPositive(expr):
 302         return sin(expr)
 303     else:
 304         return -sin(-expr)
 305 
 306 
 307 def positiveCos(expr):
 308     if pseudoPositive(expr):
 309         return cos(expr)
 310     else:
 311         return cos(-expr)
 312 
 313 
 314 def positiveTan(expr):
 315     if pseudoPositive(expr):
 316         return tan(expr)
 317     else:
 318         return -tan(-expr)
 319 
 320 
 321 def positiveCot(expr):
 322     if pseudoPositive(expr):
 323         return cot(expr)
 324     else:
 325         return -cot(-expr)
 326 
 327 
 328 def positiveCsc(expr):
 329     if pseudoPositive(expr):
 330         return csc(expr)
 331     else:
 332         return -csc(-expr)
 333 
 334 
 335 def positiveSec(expr):
 336     if pseudoPositive(expr):
 337         return sec(expr)
 338     else:
 339         return sec(-expr)
 340 
 341 
 342 def positiveSinh(expr):
 343     if pseudoPositive(expr):
 344         return sinh(expr)
 345     else:
 346         return -sinh(-expr)
 347 
 348 
 349 def positiveCosh(expr):
 350     if pseudoPositive(expr):
 351         return cosh(expr)
 352     else:
 353         return cosh(-expr)
 354 
 355 
 356 def positiveTanh(expr):
 357     if pseudoPositive(expr):
 358         return tanh(expr)
 359     else:
 360         return -tanh(-expr)
 361 
 362 
 363 def positiveCoth(expr):
 364     if pseudoPositive(expr):
 365         return coth(expr)
 366     else:
 367         return -coth(-expr)
 368 
 369 
 370 def positiveCsch(expr):
 371     if pseudoPositive(expr):
 372         return csch(expr)
 373     else:
 374         return -csch(-expr)
 375 
 376 
 377 def positiveSech(expr):
 378     if pseudoPositive(expr):
 379         return sech(expr)
 380     else:
 381         return sech(-expr)
 382 
 383 
 384 def positiveArcsin(expr):
 385     if pseudoPositive(expr):
 386         return arcsin(expr)
 387     else:
 388         return -arcsin(-expr)
 389 
 390 
 391 def positiveArctan(expr):
 392     if pseudoPositive(expr):
 393         return arctan(expr)
 394     else:
 395         return -arctan(-expr)
 396 
 397 
 398 def positiveArccot(expr):
 399     if pseudoPositive(expr):
 400         return arccot(expr)
 401     else:
 402         return -arccot(-expr)
 403 
 404 
 405 def positiveArccsc(expr):
 406     if pseudoPositive(expr):
 407         return arccsc(expr)
 408     else:
 409         return -arccsc(-expr)
 410 
 411 
 412 def positiveArcsinh(expr):
 413     if pseudoPositive(expr):
 414         return arcsinh(expr)
 415     else:
 416         return -arcsinh(-expr)
 417 
 418 
 419 def positiveArctanh(expr):
 420     if pseudoPositive(expr):
 421         return arctanh(expr)
 422     else:
 423         return -arctanh(-expr)
 424 
 425 
 426 def positiveArccoth(expr):
 427     if pseudoPositive(expr):
 428         return arccoth(expr)
 429     else:
 430         return -arccoth(-expr)
 431 
 432 
 433 def positiveArccsch(expr):
 434     if pseudoPositive(expr):
 435         return arccsch(expr)
 436     else:
 437         return -arccsch(-expr)
 438 
 439 
 440 def exp2sinhcosh(opor, opands, filter):
 441     """
 442     exp(x) => sinh(x) + cosh(x)
 443     """
 444     if opor == exp:
 445         if not filter(opands[0]):
 446             return opor(*opands)
 447         if pseudoPositive(opands[0]):
 448             return (cosh(opands[0]) + sinh(opands[0]))
 449         else:
 450             return (cosh(-opands[0]) - sinh(-opands[0]))
 451     return opor(*opands)
 452 
 453 
 454 def exp2sincos(opor, opands, filter):
 455     """
 456     exp(x) => sinh(x) + cosh(x)
 457     """
 458     if opor == exp:
 459         if not filter(opands[0]):
 460             return opor(*opands)
 461         co = I * opands[0]
 462         if pseudoPositive(co):
 463             return cos(co) - I * sin(co)
 464         else:
 465             return cos(-co) + I * sin(-co)
 466     return opor(*opands)
 467 
 468 
 469 def sinhcosh2exp(opor, opands, filter):
 470     """
 471     exp(x) => sinh(x) + cosh(x)
 472     """
 473     if opor == sinh:
 474         if not filter(opands[0]):
 475             return opor(*opands)
 476         return (exp(opands[0]) - exp(-opands[0])) / 2
 477     elif opor == cosh:
 478         if not filter(opands[0]):
 479             return opor(*opands)
 480         return (exp(opands[0]) + exp(-opands[0])) / 2
 481     return opor(*opands)
 482 
 483 
 484 def sincos2exp(opor, opands, filter):
 485     """
 486     exp(x) => sinh(x) + cosh(x)
 487     """
 488     if opor == sin:
 489         if not filter(opands[0]):
 490             return opor(*opands)
 491         return -I / 2 * (exp(I * opands[0]) - exp(-I * opands[0]))
 492     elif opor == cos:
 493         if not filter(opands[0]):
 494             return opor(*opands)
 495         return (exp(I * opands[0]) + exp(-I * opands[0])) / 2
 496     return opor(*opands)
 497 
 498 
 499 def trigo2sincos(opor, opands, filter):
 500     """
 501     [tan(x)|cot(x)] => [sin(x)/cos(x)|cos(x)/sin(x)]
 502     """
 503     if opor == tan:
 504         if not filter(opands[0]):
 505             return opor(*opands)
 506         return sin(opands[0]) / cos(opands[0])
 507     elif opor == cot:
 508         if not filter(opands[0]):
 509             return opor(*opands)
 510         return cos(opands[0]) / sin(opands[0])
 511     elif opor == csc:
 512         if not filter(opands[0]):
 513             return opor(*opands)
 514         return 1 / sin(opands[0])
 515     elif opor == sec:
 516         if not filter(opands[0]):
 517             return opor(*opands)
 518         return 1 / cos(opands[0])
 519     return opor(*opands)
 520 
 521 
 522 def trigh2sinhcosh(opor, opands, filter):
 523     """
 524     [tanh(x)|coth(x)] => [sinh(x)/cosh(x)|cosh(x)/sinh(x)]
 525     """
 526     if opor == tanh:
 527         if not filter(opands[0]):
 528             return opor(*opands)
 529         return sinh(opands[0]) / cosh(opands[0])
 530     elif opor == coth:
 531         if not filter(opands[0]):
 532             return opor(*opands)
 533         return cosh(opands[0]) / sinh(opands[0])
 534     elif opor == csch:
 535         if not filter(opands[0]):
 536             return opor(*opands)
 537         return 1 / sinh(opands[0])
 538     elif opor == sech:
 539         if not filter(opands[0]):
 540             return opor(*opands)
 541         return 1 / cosh(opands[0])
 542     return opor(*opands)
 543 
 544 
 545 def squareInPow(opor, opands, filter):
 546     if (opor != operator.pow):
 547         return(opor(*opands))
 548     if not filter[1](opands[0]):
 549         return opor(*opands)
 550     expo = opands[1]
 551     if not (expo._is_real()):
 552         return(opor(*opands))
 553     opo = opands[0].operator()
 554     if not opo:
 555         return(opor(*opands))
 556     if opo != filter[0]:
 557         return(opor(*opands))
 558     opa = opands[0].operands()[0]
 559     expo1 = floor(expo / 2)
 560     expo2 = expo - 2*expo1
 561     return filter[2](opa)**(expo2) * filter[3](opa)**(expo1)
 562 
 563 
 564 def exp2trig(opor, opands, filter):
 565     """
 566     exp(a+i*b) => (cosh(a)+sinh(a))*(cos(b)+i*sin(b))
 567     """
 568     if opor == exp:
 569         if not filter(opands[0]):
 570             return opor(*opands)
 571         res = pseudoRealAndImag(opands[0])
 572         if pseudoPositive(res[0]):
 573             pr = cosh(res[0]) + sinh(res[0])
 574         else:
 575             pr = cosh(-res[0]) - sinh(-res[0])
 576         if pseudoPositive(res[1]):
 577             pi = cos(res[1]) + I * sin(res[1])
 578         else:
 579             pi = cos(-res[1]) - I * sin(-res[1])
 580         return pr * pi
 581     return opor(*opands)
 582 
 583 
 584 def lessIinExp(opor, opands, filter):
 585     """
 586     exp(a+i*b) => exp(a)*(cos(b)+i*sin(b))
 587     """
 588     # duplicate of below ?
 589     if opor == exp:
 590         if not filter(opands[0]):
 591             return opor(*opands)
 592         res = pseudoRealAndImag(opands[0])
 593         return exp(res[0]) * (cos(res[1]) + I * sin(res[1]))
 594     return opor(*opands)
 595 
 596 
 597 def lessIinExp(opor, opands, filter):
 598     """
 599     exp(a+i*b) => exp(a)*(cos(b)+i*sin(b))
 600     """
 601     # duplicate of above ?
 602     if (opor == exp) and (opands[0].has(I)):
 603         if not filter(opands[0]):
 604             return opor(*opands)
 605         res = pseudoRealAndImag(opands[0])
 606         if pseudoPositive(res[1]):
 607             return exp(res[0]) * (cos(res[1]) + I * sin(res[1]))
 608         else:
 609             return exp(res[0]) * (cos(-res[1]) - I * sin(-res[1]))
 610     return (opor(*opands))
 611 
 612 
 613 def lessIinTrig(opor, opands, filter):
 614     """
 615     sin(i*x) => i*sinh(x) sinh, cos, cosh, tan, tanh, cot, coth
 616     """
 617     t1 = opands[0].has(I)
 618     if not t1:
 619         return (opor(*opands))
 620     if opor in [sin, cos, tan, cot, csc, sec,
 621                 sinh, cosh, tanh, coth, csch, sech]:
 622         if not filter(opands[0]):
 623             return opor(*opands)
 624         res = pseudoRealAndImag(opands[0])
 625         if res[0] != 0:
 626             return (opor(*opands))
 627     if opor == sin:
 628         return I * positiveSin(res[1])
 629     elif opor == cos:
 630         return positiveCosh(res[1])
 631     elif opor == tan:
 632         return I * positiveTanh(res[1])
 633     elif opor == cot:
 634         return -I * positiveCoth(res[1])
 635     elif opor == sec:
 636         return positiveSech(res[1])
 637     elif opor == csc:
 638         return -I / positiveCsch(res[1])
 639     elif opor == sinh:
 640         return I * positiveSin(res[1])
 641     elif opor == cosh:
 642         return positiveCos(res[1])
 643     elif opor == tanh:
 644         return I * positiveTan(res[1])
 645     elif opor == coth:
 646         return -I * positiveCot(res[1])
 647     elif opor == sech:
 648         return positiveSech(res[1])
 649     elif opor == csch:
 650         return -I / positiveCsc(res[1])
 651     return opor(*opands)
 652 
 653 
 654 def sin2tancos(opor, opands, filter):
 655     if opor == sin:
 656         if not filter(opands[0]):
 657             return opor(*opands)
 658         co = opands[0]
 659         if pseudoPositive(co):
 660             return cos(co) * tan(co)
 661         else:
 662             return - cos(-co) * tan(-co)
 663     return opor(*opands)
 664 
 665 
 666 def sinh2tanhcosh(opor, opands, filter):
 667     """exp(x) => sinh(x) + cosh(x)"""
 668     if opor == sinh:
 669         if not filter(opands[0]):
 670             return opor(*opands)
 671         co = opands[0]
 672         if pseudoPositive(co):
 673             return cosh(co) * tanh(co)
 674         else:
 675             return - cosh(-co) * tanh(-co)
 676     return opor(*opands)
 677 
 678 
 679 def sincos2tanHalf(opor, opands, filter):
 680     co = opands[0] / 2
 681     if opor == sin:
 682         if not filter(opands[0]):
 683             return opor(*opands)
 684         if pseudoPositive(co / 2):
 685             return 2*tan(co) / (1+tan(co)**2)
 686         else:
 687             return 2*tan(-co) / (1+tan(-co)**2)
 688     if opor == cos:
 689         if pseudoPositive(co / 2):
 690             return (1 - tan(co)**2) / (1 + tan(co)**2)
 691         else:
 692             return (1 - tan(-co)**2) / (1 + tan(-co)**2)
 693     return opor(*opands)
 694 
 695 
 696 def sinhcosh2tanhHalf(opor, opands, filter):
 697     """exp(x) => sinh(x) + cosh(x)"""
 698     co = opands[0] / 2
 699     if opor == sinh:
 700         if not filter(opands[0]):
 701             return opor(*opands)
 702         if pseudoPositive(co / 2):
 703             return 2*tanh(co) / (1-tanh(co)**2)
 704         else:
 705             return 2*tanh(-co) / (1-tanh(-co)**2)
 706     if opor == cosh:
 707         if not filter(opands[0]):
 708             return opor(*opands)
 709         if pseudoPositive(co / 2):
 710             return (1 + tanh(co)**2) / (1 - tanh(co)**2)
 711         else:
 712             return (1 + tanh(-co)**2) / (1 - tanh(-co)**2)
 713     return opor(*opands)
 714 
 715 
 716 def pseudoParity(opor, opands, filter):
 717     if opor == sin:
 718         if not filter(opands[0]):
 719             return opor(*opands)
 720         return positiveSin(opands[0])
 721     elif opor == cos:
 722         if not filter(opands[0]):
 723             return opor(*opands)
 724         return positiveCos(opands[0])
 725     elif opor == tan:
 726         if not filter(opands[0]):
 727             return opor(*opands)
 728         return positiveTan(opands[0])
 729     elif opor == cot:
 730         if not filter(opands[0]):
 731             return opor(*opands)
 732         return positiveCot(opands[0])
 733     elif opor == sec:
 734         if not filter(opands[0]):
 735             return opor(*opands)
 736         return positiveSec(opands[0])
 737     elif opor == csc:
 738         if not filter(opands[0]):
 739             return opor(*opands)
 740         return positiveCsc(opands[0])
 741 
 742     elif opor == sinh:
 743         if not filter(opands[0]):
 744             return opor(*opands)
 745         return positiveSinh(opands[0])
 746     elif opor == cosh:
 747         if not filter(opands[0]):
 748             return opor(*opands)
 749         return positiveCosh(opands[0])
 750     elif opor == tanh:
 751         if not filter(opands[0]):
 752             return opor(*opands)
 753         return positiveTanh(opands[0])
 754     elif opor == coth:
 755         if not filter(opands[0]):
 756             return opor(*opands)
 757         return positiveCoth(opands[0])
 758     elif opor == sech:
 759         if not filter(opands[0]):
 760             return opor(*opands)
 761         return positiveSech(opands[0])
 762     elif opor == csch:
 763         if not filter(opands[0]):
 764             return opor(*opands)
 765         return positiveCsch(opands[0])
 766 
 767     elif opor == arcsin:
 768         if not filter(opands[0]):
 769             return opor(*opands)
 770         return positiveArcsin(opands[0])
 771     elif opor == arctan:
 772         if not filter(opands[0]):
 773             return opor(*opands)
 774         return positiveArctan(opands[0])
 775     elif opor == arccot:
 776         if not filter(opands[0]):
 777             return opor(*opands)
 778         return positiveArccot(opands[0])
 779     elif opor == arccsc:
 780         if not filter(opands[0]):
 781             return opor(*opands)
 782         return positiveArccsc(opands[0])
 783 
 784     elif opor == arcsinh:
 785         if not filter(opands[0]):
 786             return opor(*opands)
 787         return positiveArcsinh(opands[0])
 788     elif opor == arctanh:
 789         if not filter(opands[0]):
 790             return opor(*opands)
 791         return positiveArctanh(opands[0])
 792     elif opor == arccoth:
 793         if not filter(opands[0]):
 794             return opor(*opands)
 795         return positiveArccoth(opands[0])
 796     elif opor == arccsch:
 797         if not filter(opands[0]):
 798             return opor(*opands)
 799         return positiveArccsch(opands[0])
 800 
 801     else:
 802         return opor(*opands)
 803 
 804 
 805 def pow2exp(opor, opands, filter):
 806     if opor == operator.pow:
 807         if not filter(opands[0], opands[1]):
 808             return opor(*opands)
 809         return exp(log(opands[0]) * opands[1])
 810     else:
 811         return opor(*opands)
 812 
 813 
 814 def exp2pow(opor, opands, filter):
 815     if opor == exp:
 816         if not filter(opands[0]):
 817             return opor(*opands)
 818         res = opands[0].match(SR.wild(0) * log(SR.wild(1)))
 819         if res:
 820             return res[SR.wild(1)]**res[SR.wild(0)]
 821     return opor(*opands)
 822 
 823 
 824 def arcsin2arccos(opor, opands, filter):
 825     """
 826     arcsin(x) => Pi/2-arccos(x)
 827     """
 828     if opor == arcsin:
 829         if not filter(opands[0]):
 830             return opor(*opands)
 831         return pi / 2 - arccos(opands[0])
 832     return opor(*opands)
 833 
 834 
 835 def arccos2arcsin(opor, opands, filter):
 836     """
 837     arccos(x) => Pi/2-arcsin(x)
 838     """
 839     if opor == arccos:
 840         if not filter(opands[0]):
 841             return opor(*opands)
 842         return pi / 2 - arcsin(opands[0])
 843     return opor(*opands)
 844 
 845 
 846 def arctrigh2log(opor, opands, filter):
 847     if opor == arcsinh:
 848         if not filter(opands[0]):
 849             return opor(*opands)
 850         return log(opands[0] + sqrt(opands[0]**2+1))
 851     elif opor == arccosh:
 852         if not filter(opands[0]):
 853             return opor(*opands)
 854         return log(opands[0] + sqrt(opands[0]**2-1))
 855     elif opor == arctanh:
 856         if not filter(opands[0]):
 857             return opor(*opands)
 858         return 1 / 2 * (log(1+opands[0])-log(1-opands[0]))
 859     elif opor == arccoth:
 860         if not filter(opands[0]):
 861             return opor(*opands)
 862         return 1 / 2 * (log(1+1/opands[0])-log(1 - 1/opands[0]))
 863     elif opor == arccsch:
 864         if not filter(opands[0]):
 865             return opor(*opands)
 866         return log(1/opands[0]+sqrt(1/opands[0]**2+1))
 867     elif opor == arcsech:
 868         if not filter(opands[0]):
 869             return opor(*opands)
 870         return log(1/opands[0]+sqrt(1/opands[0]**2-1))
 871     return opor(*opands)
 872 
 873 
 874 def arctrigo2log(opor, opands, filter):
 875     if opor == arcsin:
 876         if not filter(opands[0]):
 877             return opor(*opands)
 878         return -I * log(I * opands[0] + sqrt(1 - opands[0]**2))
 879     elif opor == arccos:
 880         if not filter(opands[0]):
 881             return opor(*opands)
 882         return -I * log(opands[0] + I * sqrt(1 - opands[0]**2))
 883     elif opor == arctan:
 884         if not filter(opands[0]):
 885             return opor(*opands)
 886         return I / 2 * (log(1 - I * opands[0]) - log(1 + I * opands[0]))
 887     elif opor == arccot:
 888         if not filter(opands[0]):
 889             return opor(*opands)
 890         return I / 2 * (log(1 - I / opands[0]) - (1 + I / opands[0]))
 891     elif opor == arccsc:
 892         if not filter(opands[0]):
 893             return opor(*opands)
 894         return -I * log(I / opands[0] + sqrt(1 - 1 / opands[0]**2))
 895     elif opor == arcsec:
 896         if not filter(opands[0]):
 897             return opor(*opands)
 898         return -I * log(1 / opands[0] + I * sqrt(1 - 1 / opands[0]**2))
 899     return opor(*opands)
 900 
 901 
 902 def log2arctanh(opor, opands, filter):
 903     if opor == ln:  # log fails
 904         if not filter(opands[0]):
 905             return opor(*opands)
 906         num = opands[0].numerator()
 907         den = opands[0].denominator()
 908         return 2 * positiveArctanh((num - den) / (num + den))
 909     return opor(*opands)
 910 
 911 
 912 def log2arctan(opor, opands, filter):
 913     if opor == ln:  # log fails
 914         if not filter(opands[0]):
 915             return opor(*opands)
 916         num = opands[0].numerator()
 917         den = opands[0].denominator()
 918         return -2 * I * positiveArctan(I * (num - den) / (num + den))
 919     return opor(*opands)
 920 
 921 
 922 def factorial2gamma(opor, opands, filter):
 923     """
 924     factorial(n) => Gamma(n+1)
 925     """
 926     if opor == factorial:
 927         if not filter(opands[0]):
 928             return opor(*opands)
 929         return gamma(opands[0] + 1)
 930     return opor(*opands)
 931 
 932 
 933 def gamma2factorial(opor, opands, filter):
 934     """
 935     Gamma(n) => factorial(n-1)
 936     """
 937     if opor == gamma:
 938         if not filter(opands[0]):
 939             return opor(*opands)
 940         return factorial(opands[0] - 1)
 941     return opor(*opands)
 942 
 943 
 944 def logMainBranch(opor, opands, filter):
 945     if opor != ln:
 946         return(opor(*opands))
 947     opo = opands[0].operator()
 948     opa = opands[0].operands()
 949     if opo != exp:
 950         return(opor(*opands))
 951     if not filter(opa[0]):
 952         return opor(*opands)
 953     else:
 954         return opa[0]
 955     return(opor(*opands))
 956 
 957 # arctrigBranchCut: [arccos|arcsin] o [+|-] [sin|cos] = 8 cas
 958 #                  [arctan|arccot] o [+|-] o [id|inverse] [tan|cot] = 16 cas
 959 
 960 
 961 def arctrigoMainBranch(opor, opands, filter):
 962     if opor == arccos:
 963         opo = opands[0].operator()
 964         opa = opands[0].operands()
 965         if opo == cos:
 966             if filter(opa[0]):
 967                 return opa[0]
 968         if opo == sin:
 969             if filter(opa[0]):
 970                 return pi / 2 - opa[0]
 971         opands2 = - opands[0]
 972         opo = opands2.operator()
 973         opa = opands2.operands()
 974         if opo == cos:
 975             if filter(opa[0]):
 976                 return pi - opa[0]
 977         if opo == sin:
 978             if filter(opa[0]):
 979                 return pi / 2 - opa[0]
 980         return opor(*opands)
 981 
 982     if opor == arcsin:
 983         opo = opands[0].operator()
 984         opa = opands[0].operands()
 985         if opo == sin:
 986             if filter(opa[0]):
 987                 return opa[0]
 988         if opo == cos:
 989             if filter(opa[0]):
 990                 return pi / 2 - opa[0]
 991         opands2 = - opands[0]
 992         opo = opands2.operator()
 993         opa = opands2.operands()
 994         if opo == sin:
 995             if filter(opa[0]):
 996                 return - opa[0]
 997         if opo == cos:
 998             if filter(opa[0]):
 999                 return opa[0] - pi / 2
1000         return opor(*opands)
1001 
1002     if opor == arctan:
1003         opo = opands[0].operator()
1004         opa = opands[0].operands()
1005         if opo == tan:
1006             if filter(opa[0]):
1007                 return opa[0]
1008         if opo == cot:
1009             if filter(opa[0]):
1010                 return pi / 2 - opa[0]
1011         opands2 = - opands[0]
1012         opo = opands2.operator()
1013         opa = opands2.operands()
1014         if opo == tan:
1015             if filter(opa[0]):
1016                 return - opa[0]
1017         if opo == cot:
1018             if filter(opa[0]):
1019                 return opa[0] - pi / 2
1020         opands2 = 1/opands[0]
1021         opo = opands2.operator()
1022         opa = opands2.operands()
1023         if opo == cot:
1024             if filter(opa[0]):
1025                 return opa[0]
1026         if opo == tan:
1027             if filter(opa[0]):
1028                 return pi / 2 - opa[0]
1029         opands2 = -opands2
1030         opo = opands2.operator()
1031         opa = opands2.operands()
1032         if opo == cot:
1033             if filter(opa[0]):
1034                 return -opa[0]
1035         if opo == tan:
1036             if filter(opa[0]):
1037                 return opa[2] - pi / 2
1038         return opor(*opands)
1039 
1040     if opor == arccot:
1041         opo = opands[0].operator()
1042         opa = opands[0].operands()
1043         if opo == cot:
1044             if filter(opa[0]):
1045                 return opa[0]
1046         if opo == tan:
1047             if filter(opa[0]):
1048                 return pi / 2 - opa[0]
1049         opands2 = - opands[0]
1050         opo = opands2.operator()
1051         opa = opands2.operands()
1052         if opo == cot:
1053             if filter(opa[0]):
1054                 return - opa[0]
1055         if opo == tan:
1056             if filter(opa[0]):
1057                 return opa[0] - pi / 2
1058         opands2 = 1 / opands[0]
1059         opo = opands2.operator()
1060         opa = opands2.operands()
1061         if opo == tan:
1062             if filter(opa[0]):
1063                 return opa[0]
1064         if opo == cot:
1065             if filter(opa[0]):
1066                 return pi / 2 - opa[0]
1067         opands2 = -opands2
1068         opo = opands2.operator()
1069         opa = opands2.operands()
1070         if opo == tan:
1071             if filter(opa[0]):
1072                 return -opa[0]
1073         if opo == cot:
1074             if filter(opa[0]):
1075                 return opa[2] - pi / 2
1076         return opor(*opands)
1077 
1078     return opor(*opands)