attachment:rewrite-20110123.sage of symbolics/rewrite

   1 # François Maltey - october 2010 
   2 #
   3 # The "mapexpression" function maps recursively an expression as a tree.
   4 #
   5 # This function 
   6 #    remains the framework of the expression, 
   7 #    treats first the leaves and then the composed expression, and
   8 #    can change only some subtree.
   9 # This change can be done 
  10 #    at the main level only, 
  11 #    at some levels in the tree,
  12 #    at all the levels.
  13 # 
  14 # The function goes thru sum and product and can consider 
  15 # that the depth of the subexpression doesn't change.
  16 # 
  17 # expr = the current expression
  18 # mulDepth = 0 or 1, if 0 the depth in the tree remains the same in a sum.
  19 # addDepth = 0 or 1, if 1 the depth in the tree increases in a product.
  20 # fct = the effective function called for subtrees 
  21 # param = the parameter of fct
  22 # level = -1 for a fully recursive call, or the list of levels to treat.
  23 # level = [0] for the first main level.
  24 
  25 import operator
  26 
  27 def mapexpression (expr, fct, param, level, addDepth=0, mulDepth=0) :
  28   def mapex (expr, depth) :               # a very local function
  29     opor = expr.operator()     
  30     opands = expr.operands()  
  31     if (opor == None) : return expr       # a leaf in the expression tree 
  32     if (opor == operator.add) :           # recursive call thru sum
  33       opands = map (lambda ex : mapex (ex, depth+addDepth), opands)
  34       return add (opands)
  35     if (opor == operator.mul) :           # recursive call thru mul
  36       opands = map (lambda ex : mapex (ex, depth+mulDepth), opands)
  37       return mul (opands)   
  38     if (level == -1) or (level[-1] >= depth) : # recursive call over operands
  39       opands = map (lambda ex : mapex (ex, depth+1), opands)
  40     if level == -1 or depth in level :    # root of the subtree must be changed
  41       return fct (opor, opands, param)
  42     return opor (*opands) # opands may or maynot be changed by a recursive call
  43   return mapex (expr, 0) 
  44 
  45 #
  46 # Choose the sign of expressions in respect of the opposite form.
  47 # If a-b is said "positive" then -(a-b) is negative.
  48 # For 
  49 
  50 def pseudoPositive (expr) :
  51   if expr._is_real() : 
  52     return bool (RR(expr) >= 0) # can be improved by ._is_positive()
  53   if expr._is_numeric () :      # see ._is_positive()
  54     return bool ((expr.real() > 0) or (expr.real() == 0 and expr.imag() > 0))
  55   if expr._is_symbol() : return True  # a_variable as x or a "is positive"
  56   opor = expr.operator()
  57   opands = expr.operands()
  58   if opor == operator.mul :            # look at the last factor in a product
  59     return pseudoPositive (opands[-1]) # it's the number if there is one
  60   if opor == operator.add :    # look at the first term in a sum
  61     return pseudoPositive (opands[0]) 
  62   return True    # for functions call as sin(x)...
  63 
  64 def pseudoRealAndImag (expr) : # try to get (a,b) from a+i*b
  65   opands = expr.operands()     # don't decompose (a+i*b)*(x+i*y)
  66   opor = expr.operator()       
  67   if opor == operator.mul :    # but treat Complex_number * expressions   
  68     coeff = 1 
  69     unit = 1
  70     for ex in opands :
  71       if ex._is_numeric() : coeff = coeff * ex
  72       else : unit = unit*ex
  73     return (coeff.real()*unit, coeff.imag()*unit)
  74   elif opor == operator.add :  # and treat sum
  75     rp = 0
  76     ip = 0
  77     for ex in opands :
  78       if ex._is_numeric() : rp = rp + ex.real() ; ip = ip + ex.imag()
  79       res = pseudoRealAndImag(ex)
  80       rp = rp + res[0] 
  81       ip = ip + res[1]
  82     return (rp, ip)
  83   return (expr, 0)             # in doubt remain the expression as real
  84 
  85 # There are 2 uses of rewrite : 
  86 #   either by the 2 parameters source=... and target=...
  87 #   either by the name "source2target" of this rule
  88 # This function transforms the first method to the second one
  89 
  90 def searchRule (source, target) :
  91   if source in [sin,cos] and target in [sinh,cosh] : return ["trigo2trigh"]
  92   elif source == exp and target in [sin,cos] : return ["exp2trigo"]
  93   elif source == exp and target in [sinh,cosh] : return ["exp2trigh"]
  94   elif source in [sin,cos] and target == exp : return ["exp2trigo"]
  95   elif source in [sinh,cosh] and target == exp : return ["exp2trigh"]
  96   else : return "can't match source and target, use one rule or a list of rules"
  97 
  98 def rewrite (expr, rule=None, source=None, target=None, filter=None, level=-1) :
  99   if rule==target==None :    # the function controls the most common errors 
 100     return "must choose either rule=... or target=..."
 101   if rule!=None and target!=None :
 102     return "must choose either rule=... or target=..."
 103   if target==None and source!= None :
 104     return "must choose target=... for a definded source"
 105   if rule==None : rule=searchRules (source, target)
 106   elif type(rule)==str : rule=[rule]
 107   else : return 'a rule must be a string or a list of strings'
 108   if filter==None : filter = lambda ex : true
 109   if type(filter).__name__ != 'function' :
 110     return 'filter must be a function'
 111   if level != -1 and type(level) != list :
 112     return 'level must be -1 for a fully recursive call or a list of integer'
 113   k=0 
 114   while k<len(rule) : 
 115     if rule[k] == "sinhcosh2exp" :   # sinh(x) -> (exp(x)-exp(-x))/2, ...
 116       expr = mapexpression (expr, sinhcosh2exp, filter, level)
 117     elif rule[k] == "sincos2exp" :   # sin(x) -> (exp(ix)-exp(-ix))/(2i), ...
 118       expr = mapexpression (expr, sincos2exp, filter, level)
 119     elif rule[k] == "exp2sinhcosh" : # exp(x) -> sinh(x)+cosh(x)
 120       expr = mapexpression (expr, exp2sinhcosh, filter, level)
 121     elif rule[k] == "exp2sincos" :   # exp(x)=exp(i*(-ix)) -> cos(ix)-i*sin(ix) 
 122       expr = mapexpression (expr, exp2sincos, filter, level)
 123 
 124     elif rule[k] == "trigo2trigh" :  # cos(x) -> cosh(ix), ...
 125       expr = mapexpression (expr, trigo2trigh, filter, level)
 126     elif rule[k] == "trigh2trigo" :  # cosh(x) -> cos(i*x), ...
 127       expr = mapexpression (expr, trigh2trigo, filter, level)
 128 
 129     elif rule[k]== "trigo2sincos" :  # tan, cot, sec, csc -> sin, cos
 130       expr = mapexpression (expr, trigo2sincos, filter, level)
 131     elif rule[k]== "trigh2sinhcosh" :# tanh, coth, sech, csch -> sinh, cosh
 132       expr = mapexpression (expr, trigh2sinhcosh, filter, level)
 133 
 134     elif rule[k]== "exp2trig" : #exp(a+i*b)->(cosh(a)+sinh(a))*(cos(b)+i*sin(b))
 135       expr = mapexpression (expr, exp2trig, filter, level)
 136     elif rule[k]== "lessIinExp" :    # exp(a+i*b) -> exp(a)*(cos(b)+i*sin(b))
 137       expr = mapexpression (expr, lessIinExp, filter, level)
 138     elif rule[k]== "lessIinTrig" :   # cos(i*x) -> cosh(x), ...
 139       expr = mapexpression (expr, lessIinTrig,filter,level)
 140     
 141     elif rule[k]== "trigo2exp" : rule.extend(["sincos2exp", "trigo2sincos"])
 142     elif rule[k]== "trigh2exp" : rule.extend(["sinhcosh2exp", "trigh2sinhcosh"])
 143     elif rule[k]== "trig2exp" : rule.extend(["trigo2exp", "trigh2exp"])
 144     
 145     elif rule[k]== "cos22sin" :
 146       expr = \
 147         mapexpression (expr, squareInPow,[cos, filter, positiveCos, \
 148                                       lambda ex:1-positiveSin(ex)^2], level)
 149     elif rule[k]== "sin22cos" :
 150       expr = \
 151         mapexpression (expr, \
 152           squareInPow, [sin, filter, positiveSin, lambda ex:1-positiveCos(ex)^2], level)
 153     elif rule[k]== "cosh22sinh" :
 154       expr = \
 155         mapexpression (expr, \
 156           squareInPow, [cosh, filter, positiveCosh, lambda ex:1+positiveSinh(ex)^2], level)
 157     elif rule[k]== "sinh22cosh" :
 158       expr = \
 159         mapexpression (expr, \
 160           squareInPow, [sinh, filter, positiveSinh, lambda ex:positiveCosh(ex)^2-1], level)
 161     elif rule[k]== "tan22cos" :
 162       expr = \
 163         mapexpression (expr,\
 164           squareInPow, [tan, filter, positiveTan, lambda ex:1/positiveCos(ex)^2-1], level)
 165     elif rule[k]== "cot22sin" :
 166       expr = \
 167         mapexpression (expr,\
 168           squareInPow, [cot, filter, positiveCot, lambda ex:1/positiveSin(ex)^2-1], level)
 169     elif rule[k]== "tanh22cosh" :
 170       expr = \
 171         mapexpression (expr, \
 172           squareInPow, [tanh,filter,positiveTanh,lambda ex:1-1/positiveCosh(ex)^2], level)
 173     elif rule[k]== "coth22sinh" :
 174       expr = \
 175         mapexpression (expr, \
 176           squareInPow, [coth,filter,positiveCoth,lambda ex:1/positiveSinh(ex)^2-1], level)
 177     elif rule[k]== "cos22tan" :
 178       expr = mapexpression (expr, squareInPow, \
 179        [cos, filter, positiveCos, lambda ex:1/(positiveTan(ex)^2+1)], level)
 180     elif rule[k]== "tancot22sincos" : rule.extend(["tan22cos", "cot22sin"])
 181     elif rule[k]== "tanhcoth22sinhcosh" : rule.extend(["tanh22cosh", "coth22sinh"])
 182 
 183     elif rule[k]== "sin2tancos" :
 184       expr = mapexpression (expr, sin2tancos, filter, level)
 185     elif rule[k] == "sincos2tan" :
 186       expr = mapexpression (expr, sin2tancos, filter, level)
 187       expr = expr.rational_simplify()
 188       expr = mapexpression (expr, squareInPow, \
 189        [cos, filter, positiveCos, lambda ex:1/(positiveTan(ex)^2+1)], level)
 190     elif rule[k]== "sinh2tanhcosh" :
 191       expr = mapexpression (expr, sinh2tanhcosh, filter, level)
 192     
 193     elif rule[k]== "sincos2tanHalf" :
 194       expr = mapexpression (expr, sincos2tanHalf, filter, level)
 195     elif rule[k]== "sinhcosh2tanhHalf" :
 196       expr = mapexpression (expr, sinhcosh2tanhHalf, filter, level)
 197 
 198     elif rule[k]== "exp2pow" :
 199       expr = mapexpression (expr, exp2pow, filter, level)
 200     elif rule[k]== "pow2exp" :
 201       expr = mapexpression (expr, pow2exp, filter, level)
 202 
 203 
 204     elif rule[k]== "arcsin2arccos" :
 205       expr = mapexpression (expr, arcsin2arccos, filter, level)
 206     elif rule[k]== "arccos2arcsin" :
 207       expr = mapexpression (expr, arccos2arcsin, filter, level)
 208     elif rule[k]== "arctrigo2log" :
 209       expr = mapexpression (expr, arctrigo2log, filter, level)
 210     elif rule[k]== "arctrigh2log" :
 211       expr = mapexpression (expr, arctrigh2log, filter, level)
 212     elif rule[k]== "log2arctan" :
 213       expr = mapexpression (expr, log2arctan, filter, level)
 214     elif rule[k]== "log2arctanh" :
 215       expr = mapexpression (expr, log2arctanh, filter, level)
 216     elif rule[k]== "factorial2gamma" :
 217       expr = mapexpression (expr, factorial2gamma, filter, level)
 218     elif rule[k]== "gamma2factorial" :
 219       expr = mapexpression (expr, gamma2factorial, filter, level)
 220 #    elif rule[k]== "binomial2factorial" :
 221 #      expr = mapexpression (expr, binomial2factorial, filter, level)
 222 
 223     elif rule[k]== "pseudoParity" :
 224       expr = mapexpression (expr, pseudoParity, filter, level)
 225 
 226     elif rule[k]== "logMainBranch" :
 227       expr = mapexpression (expr, logMainBranch, filter, level)
 228     elif rule[k]== "arctrigoMainBranch" :
 229       expr = mapexpression (expr, arctrigoMainBranch, filter, level)
 230 
 231 
 232     else : return "unknown rule"
 233     k=k+1
 234   return expr
 235 
 236 def positiveSin (expr) : 
 237   if pseudoPositive (expr) : return sin (expr)
 238   else : return -sin(-expr)
 239 def positiveCos (expr) : 
 240   if pseudoPositive (expr) : return cos (expr)
 241   else : return cos(-expr)
 242 
 243 def positiveTan (expr) : 
 244   if pseudoPositive (expr) : return tan (expr)
 245   else : return -tan(-expr)
 246 def positiveCot (expr) : 
 247   if pseudoPositive (expr) : return cot (expr)
 248   else : return -cot(-expr)
 249 
 250 def positiveCsc (expr) : 
 251   if pseudoPositive (expr) : return csc (expr)
 252   else : return -csc(-expr)
 253 def positiveSec (expr) : 
 254   if pseudoPositive (expr) : return sec (expr)
 255   else : return sec(-expr)
 256 
 257 def positiveSinh (expr) : 
 258   if pseudoPositive (expr) : return sinh (expr)
 259   else : return -sinh(-expr)
 260 def positiveCosh (expr) : 
 261   if pseudoPositive (expr) : return cosh (expr)
 262   else : return cosh(-expr)
 263 
 264 def positiveTanh (expr) : 
 265   if pseudoPositive (expr) : return tanh (expr)
 266   else : return -tanh(-expr)
 267 def positiveCoth (expr) : 
 268   if pseudoPositive (expr) : return coth (expr)
 269   else : return -coth(-expr)
 270 
 271 def positiveCsch (expr) : 
 272   if pseudoPositive (expr) : return csch (expr)
 273   else : return -csch(-expr)
 274 def positiveSech (expr) : 
 275   if pseudoPositive (expr) : return sech (expr)
 276   else : return sech(-expr)
 277 
 278 def positiveArcsin (expr) : 
 279   if pseudoPositive (expr) : return arcsin (expr)
 280   else : return -arcsin(-expr)
 281 def positiveArctan (expr) : 
 282   if pseudoPositive (expr) : return arctan (expr)
 283   else : return -arctan(-expr)
 284 def positiveArccot (expr) : 
 285   if pseudoPositive (expr) : return arccot (expr)
 286   else : return -arccot(-expr)
 287 def positiveArccsc (expr) : 
 288   if pseudoPositive (expr) : return arccsc (expr)
 289   else : return -arccsc(-expr)
 290 def positiveArcsinh (expr) : 
 291   if pseudoPositive (expr) : return arcsinh (expr)
 292   else : return -arcsinh(-expr)
 293 def positiveArctanh (expr) : 
 294   if pseudoPositive (expr) : return arctanh (expr)
 295   else : return -arctanh(-expr)
 296 def positiveArccoth (expr) : 
 297   if pseudoPositive (expr) : return arccoth (expr)
 298   else : return -arccoth(-expr)
 299 def positiveArccsch (expr) : 
 300   if pseudoPositive (expr) : return arccsch (expr)
 301   else : return -arccsch(-expr)
 302 
 303 def exp2sinhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 304   if opor == exp :
 305     if not filter(opands[0]) : return opor (*opands)
 306     if pseudoPositive (opands[0]) :
 307       return (cosh (opands[0]) + sinh (opands[0]))
 308     else : 
 309        return (cosh (-opands[0]) - sinh (-opands[0]))
 310   return opor (*opands)
 311 
 312 def exp2sincos (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 313   if opor == exp :
 314     if not filter(opands[0]) : return opor (*opands)
 315     co = I*opands[0] 
 316     if pseudoPositive (co) : return cos(co) - I * sin(co)
 317     else : return cos(-co) + I * sin(-co)
 318   return opor (*opands)
 319 
 320 def sinhcosh2exp (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 321   if opor == sinh :
 322     if not filter(opands[0]) : return opor (*opands)
 323     return (exp (opands[0]) - exp (-opands[0]))/2
 324   elif opor == cosh :
 325     if not filter(opands[0]) : return opor (*opands)
 326     return (exp (opands[0]) + exp (-opands[0]))/2
 327   return opor (*opands)
 328 
 329 def sincos2exp (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 330   if opor == sin :
 331     if not filter(opands[0]) : return opor (*opands)
 332     return -I/2 * (exp (I*opands[0]) - exp (-I*opands[0]))
 333   elif opor == cos :
 334     if not filter(opands[0]) : return opor (*opands)
 335     return (exp (I*opands[0]) + exp (-I*opands[0]))/2
 336   return opor (*opands)
 337 
 338 def trigo2sincos (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 339   if opor == tan :
 340     if not filter(opands[0]) : return opor (*opands)
 341     return sin (opands[0]) / cos (opands[0])
 342   elif opor == cot :
 343     if not filter(opands[0]) : return opor (*opands)
 344     return cos (opands[0]) / sin (opands[0])
 345   elif opor == csc :
 346     if not filter(opands[0]) : return opor (*opands)
 347     return 1 / sin (opands[0])
 348   elif opor == sec : 
 349     if not filter(opands[0]) : return opor (*opands)
 350     return 1 / cos (opands[0])
 351   return opor (*opands)
 352 
 353 def trigh2sinhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 354   if opor == tanh :
 355     if not filter(opands[0]) : return opor (*opands)
 356     return sinh (opands[0]) / cosh (opands[0])
 357   elif opor == coth :
 358     if not filter(opands[0]) : return opor (*opands)
 359     return cosh (opands[0]) / sinh (opands[0])
 360   elif opor == csch :
 361     if not filter(opands[0]) : return opor (*opands)
 362     return 1 / sinh (opands[0])
 363   elif opor == sech :
 364     if not filter(opands[0]) : return opor (*opands)
 365     return 1 / cosh (opands[0])
 366   return opor (*opands)
 367 
 368 def squareInPow(opor, opands, filter) :
 369   if (opor != operator.pow) : return(opor (*opands))
 370   if not filter[1](opands[0]) : return opor (*opands)
 371   expo = opands[1]
 372   if not (expo._is_real()) : return(opor (*opands))
 373   opo = opands[0].operator()
 374   if not opo : return(opor (*opands))
 375   if opo != filter[0] : return(opor (*opands))
 376   opa = opands[0].operands()[0]
 377   expo1 = floor(expo/2)
 378   expo2 = expo - 2*expo1
 379   return filter[2](opa)^(expo2) * filter[3](opa)^(expo1)
 380 
 381 def lessIinExp (opor, opands, filter) :
 382   if opor == exp :
 383     if not filter(opands[0]) : return opor (*opands)
 384     res = pseudoRealAndImag (opands[0]) 
 385     return exp(res[0]) * (cos(res[1]) + I*sin(res[1]))
 386   return opor (*opands)
 387 
 388 def exp2trig (opor, opands, filter) :
 389   if opor == exp :
 390     if not filter(opands[0]) : return opor (*opands)
 391     res = pseudoRealAndImag (opands[0])
 392     if pseudoPositive (res[0]) : pr = cosh(res[0])+sinh(res[0])
 393     else : pr = cosh(-res[0])-sinh(-res[0])
 394     if pseudoPositive (res[1]) : pi = cos(res[1])+I*sin(res[1]) 
 395     else : pi = cos(-res[1])-I*sin(-res[1])
 396     return  pr * pi
 397   return opor (*opands)
 398 
 399 def lessIinExp (opor, opands, filter) :
 400   if (opor == exp) and (opands[0].has(i))  :
 401     if not filter(opands[0]) : return opor (*opands)
 402     res = pseudoRealAndImag (opands[0])
 403     if pseudoPositive (res[1]) :
 404       return exp(res[0]) * (cos(res[1]) + I*sin(res[1]))
 405     else :
 406       return exp(res[0]) * (cos(-res[1]) - I*sin(-res[1]))
 407   return (opor(*opands))
 408 
 409 def lessIinTrig (opor, opands, filter) :
 410   t1 = opands[0].has(i)
 411   if not t1 :   return (opor(*opands))
 412   if opor in [sin,cos,tan,cot,csc,sec,sinh,cosh,tanh,coth,csch,sech] :
 413     if not filter(opands[0]) : return opor (*opands)
 414     res = pseudoRealAndImag (opands[0]) 
 415     if res[0] != 0 : return (opor(*opands))
 416   if opor == sin : return i*positiveSin(res[1])
 417   elif opor == cos : return positiveCosh(res[1])
 418   elif opor == tan : return i*positiveTanh(res[1])
 419   elif opor == cot : return -i*positiveCoth(res[1])
 420   elif opor == sec : return positiveSech(res[1])
 421   elif opor == csc : return -i/positiveCsch(res[1])
 422   elif opor == sinh : return i*positiveSin(res[1])
 423   elif opor == cosh : return positiveCos(res[1])
 424   elif opor == tanh : return i*positiveTan(res[1])
 425   elif opor == coth : return -i*positiveCot(res[1])
 426   elif opor == sech : return positiveSech(res[1])
 427   elif opor == csch : return -i/positiveCsc(res[1])
 428   return (opor(*opands))
 429 
 430 def sin2tancos (opor, opands, filter)  : 
 431   if opor == sin :
 432     if not filter(opands[0]) : return opor (*opands)
 433     co = opands[0] 
 434     if pseudoPositive (co) : return cos(co) * tan(co)
 435     else : return - cos(-co) * tan(-co)
 436   return opor (*opands)
 437 
 438 def sinh2tanhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 439   if opor == sinh :
 440     if not filter(opands[0]) : return opor (*opands)
 441     co = opands[0] 
 442     if pseudoPositive (co) : return cosh(co) * tanh(co)
 443     else : return - cosh(-co) * tanh(-co)
 444   return opor (*opands)
 445 
 446 def sincos2tanHalf (opor, opands, filter)  : 
 447   co = opands[0]/2 
 448   if opor == sin :
 449     if not filter(opands[0]) : return opor (*opands)
 450     if pseudoPositive (co/2) : return 2*tan(co)/(1+tan(co)^2)
 451     else : return 2*tan(-co)/(1+tan(-c)^2)
 452   if opor == cos :
 453     if pseudoPositive (co/2) : return (1-tan(co)^2)/(1+tan(co)^2)
 454     else : return (1-tan(-co)^2)/(1+tan(-c)^2)
 455   return opor (*opands)
 456 
 457 def sinhcosh2tanhHalf (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 458   co = opands[0]/2 
 459   if opor == sinh :
 460     if not filter(opands[0]) : return opor (*opands)
 461     if pseudoPositive (co/2) : return 2*tanh(co)/(1-tanh(co)^2)
 462     else : return 2*tanh(-co)/(1-tanh(-c)^2)
 463   if opor == cosh :
 464     if not filter(opands[0]) : return opor (*opands)
 465     if pseudoPositive (co/2) : return (1+tanh(co)^2)/(1-tanh(co)^2)
 466     else : return (1+tanh(-co)^2)/(1-tanh(-c)^2)
 467   return opor (*opands)
 468 
 469 def pseudoParity (opor, opands, filter) :
 470   if opor == sin : 
 471     if not filter(opands[0]) : return opor (*opands)
 472     return positiveSin (opands[0])
 473   elif opor == cos :
 474     if not filter(opands[0]) : return opor (*opands)
 475     return positiveCos (opands[0])
 476   elif opor == tan :
 477     if not filter(opands[0]) : return opor (*opands)
 478     return positiveTan (opands[0])  
 479   elif opor == cot :
 480     if not filter(opands[0]) : return opor (*opands)
 481     return positiveCot (opands[0])  
 482   elif opor == sec :
 483     if not filter(opands[0]) : return opor (*opands)
 484     return positiveSec (opands[0])  
 485   elif opor == csc :
 486     if not filter(opands[0]) : return opor (*opands)
 487     return positiveCsc (opands[0])
 488 
 489   elif opor == sinh :
 490     if not filter(opands[0]) : return opor (*opands)
 491     return positiveSinh (opands[0])
 492   elif opor == cosh :
 493     if not filter(opands[0]) : return opor (*opands)
 494     return positiveCosh (opands[0])
 495   elif opor == tanh :
 496     if not filter(opands[0]) : return opor (*opands)
 497     return positiveTanh (opands[0])  
 498   elif opor == coth :
 499     if not filter(opands[0]) : return opor (*opands)
 500     return positiveCoth (opands[0])  
 501   elif opor == sech :
 502     if not filter(opands[0]) : return opor (*opands)
 503     return positiveSech (opands[0])  
 504   elif opor == csch :
 505     if not filter(opands[0]) : return opor (*opands)
 506     return positiveCsch (opands[0])
 507 
 508   elif opor == arcsin :
 509     if not filter(opands[0]) : return opor (*opands)
 510     return positiveArcsin (opands[0])
 511   elif opor == arctan :
 512     if not filter(opands[0]) : return opor (*opands)
 513     return positiveArctan (opands[0])  
 514   elif opor == arccot :
 515     if not filter(opands[0]) : return opor (*opands)
 516     return positiveArccot (opands[0])  
 517   elif opor == arccsc :
 518     if not filter(opands[0]) : return opor (*opands)
 519     return positiveArccsc (opands[0])  
 520 
 521   elif opor == arcsinh :
 522     if not filter(opands[0]) : return opor (*opands)
 523     return positiveArcsinh (opands[0])
 524   elif opor == arctanh :
 525     if not filter(opands[0]) : return opor (*opands)
 526     return positiveArctanh (opands[0])  
 527   elif opor == arccoth :
 528     if not filter(opands[0]) : return opor (*opands)
 529     return positiveArccoth (opands[0])  
 530   elif opor == arccsch :
 531     if not filter(opands[0]) : return opor (*opands)
 532     return positiveArccsch (opands[0])  
 533 
 534   else : return opor(*opands)
 535 
 536 def pow2exp (opor, opands, filter) :
 537   if opor == operator.pow : 
 538     if not filter(opands[0], opands[1]) : return opor (*opands)
 539     return (exp (log(opands[0]) * opands[1]))
 540   else : return opor(*opands)
 541 
 542 def exp2pow (opor, opands, filter) :
 543   if opor == exp : 
 544     if not filter(opands[0]) : return opor (*opands)
 545     res = opands[0].match(SR.wild(0)*log(SR.wild(1)))
 546     if res : return res[SR.wild(1)]^res[SR.wild(0)]
 547   return opor(*opands)
 548 
 549 def arcsin2arccos (opor, opands, filter) :
 550   if opor == arcsin : 
 551     if not filter(opands[0]) : return opor (*opands)
 552     return pi/2-arccos(opands[0])
 553   return opor(*opands)
 554 
 555 def arccos2arcsin (opor, opands, filter) :
 556   if opor == arccos : 
 557     if not filter(opands[0]) : return opor (*opands)
 558     return pi/2-arcsin(opands[0])
 559   return opor(*opands)
 560 
 561 def arctrigh2log (opor, opands, filter) :
 562   if opor == arcsinh : 
 563     if not filter(opands[0]) : return opor (*opands)
 564     return log(opands[0]+sqrt(opands[0]^2+1))
 565   elif opor == arccosh : 
 566     if not filter(opands[0]) : return opor (*opands)
 567     return log(opands[0]+sqrt(opands[0]^2-1))
 568   elif opor == arctanh :
 569     if not filter(opands[0]) : return opor (*opands)
 570     return 1/2*(log(1+opands[0])-log(1-opands[0]))
 571   elif opor == arccoth :
 572     if not filter(opands[0]) : return opor (*opands)
 573     return 1/2*(log(1+1/opands[0])-log(1-1/opands[0]))
 574   elif opor == arccsch : 
 575     if not filter(opands[0]) : return opor (*opands)
 576     return log(1/opands[0]+sqrt(1/opands[0]^2+1)) 
 577   elif opor == arcsech : 
 578     if not filter(opands[0]) : return opor (*opands)
 579     return log(1/opands[0]+sqrt(1/opands[0]^2-1))
 580   return opor(*opands)
 581 
 582 def arctrigo2log (opor, opands, filter) :
 583   if opor == arcsin : 
 584     if not filter(opands[0]) : return opor (*opands)
 585     return -i*log(i*opands[0]+sqrt(1-opands[0]^2))
 586   elif opor == arccos : 
 587     if not filter(opands[0]) : return opor (*opands)
 588     return -i*log(opands[0]+i*sqrt(1-opands[0]^2))
 589   elif opor == arctan : 
 590     if not filter(opands[0]) : return opor (*opands)
 591     return i/2*(log(1-i*opands[0])-log(1+i*opands[0]))
 592   elif opor == arccot : 
 593     if not filter(opands[0]) : return opor (*opands)
 594     return i/2*(log(1-i/opands[0])-(1+i/opands[0]))
 595   elif opor == arccsc : 
 596     if not filter(opands[0]) : return opor (*opands)
 597     return -i*log(i/opands[0]+sqrt(1-1/opands[0]^2))
 598   elif opor == arcsec : 
 599     if not filter(opands[0]) : return opor (*opands)
 600     return -i*log(1/opands[0]+i*sqrt(1-1/opands[0]^2))
 601   return opor(*opands)
 602 
 603 def log2arctanh (opor, opands, filter) :
 604   if opor == ln :  # log fails
 605     if not filter(opands[0]) : return opor (*opands)
 606     num = opands[0].numerator()
 607     den = opands[0].denominator()
 608     return 2*positiveArctanh ((num-den)/(num+den)) 
 609   return opor(*opands)
 610 
 611 def log2arctan (opor, opands, filter) :
 612   if opor == ln :  # log fails
 613     if not filter(opands[0]) : return opor (*opands)
 614     num = opands[0].numerator()
 615     den = opands[0].denominator()
 616     return -2*i*positiveArctan (i*(num-den)/(num+den)) 
 617   return opor(*opands)
 618 
 619 def factorial2gamma (opor, opands, filter) :
 620   if opor == factorial : 
 621     if not filter(opands[0]) : return opor (*opands)
 622     return gamma(opands[0]+1)
 623   return opor(*opands)
 624 def gamma2factorial (opor, opands, filter) :
 625   if opor == gamma : 
 626     if not filter(opands[0]) : return opor (*opands)
 627     return factorial(opands[0]-1)
 628   return opor(*opands)
 629 
 630 #def binomial2factorial (opor, opands, filter) :
 631 #  if opor == binomial : 
 632 #    return    factorial (opand[0]) \
 633 #           / (factorial (operands[1]) * factorial (operands[0]-operands[1]))
 634 #  return opor(*opands)
 635 
 636 def logMainBranch (opor, opands, filter) :
 637 #  if opor != log : return(opor (*opands))
 638   if opor != ln : return(opor (*opands))
 639   opo = opands[0].operator()
 640   opa = opands[0].operands()
 641   if opo != exp : return(opor (*opands))
 642   if not filter(opa[0]) : return opor (*opands)
 643   else : return opa[0]
 644   return(opor (*opands))
 645 
 646 # arctrigBranchCut : [arccos|arcsin] o [+|-] [sin|cos] = 8 cas
 647 #                  [arctan|arccot] o [+|-] o [id|inverse] [tan|cot] = 16 cas
 648 
 649 def arctrigoMainBranch (opor, opands, filter) :
 650   if opor == arccos :
 651     opo = opands[0].operator()
 652     opa = opands[0].operands()
 653     if opo == cos : 
 654       if filter(opa[0]) : return opa[0]
 655     if opo == sin : 
 656       if filter(opa[0]) : return pi/2-opa[0] 
 657     opands2 = - opands[0]
 658     opo = opands2.operator()
 659     opa = opands2.operands()
 660     if opo == cos : 
 661       if filter(opa[0]) : return pi - opa[0]
 662     if opo == sin : 
 663       if filter(opa[0]) : return pi/2 - opa[0]
 664     return opor (*opands)
 665 
 666   if opor == arcsin :
 667     opo = opands[0].operator()
 668     opa = opands[0].operands()
 669     if opo == sin : 
 670       if filter(opa[0]) : return opa[0]
 671     if opo == cos : 
 672       if filter(opa[0]) : return pi/2 - opa[0] 
 673     opands2 = - opands[0]
 674     opo = opands2.operator()
 675     opa = opands2.operands()
 676     if opo == sin : 
 677       if filter(opa[0]) : return - opa[0]
 678     if opo == cos : 
 679       if filter(opa[0]) : return opa[0] - pi/2
 680     return opor (*opands)
 681 
 682   if opor == arctan :
 683     opo = opands[0].operator()
 684     opa = opands[0].operands()
 685     if opo == tan : 
 686       if filter(opa[0]) : return opa[0]
 687     if opo == cot : 
 688       if filter(opa[0]) : return pi/2-opa[0] 
 689     opands2 = - opands[0]
 690     opo = opands2.operator()
 691     opa = opands2.operands()
 692     if opo == tan : 
 693       if filter(opa[0]) : return - opa[0]
 694     if opo == cot : 
 695       if filter(opa[0]) : return opa[0] - pi/2 
 696     opands2 = 1/opands[0]
 697     opo = opands2.operator()
 698     opa = opands2.operands()
 699     if opo == cot : 
 700       if filter(opa[0]) : return  opa[0]
 701     if opo == tan : 
 702       if filter(opa[0]) : return pi/2-opa[0] 
 703     opands2 = -opands2
 704     opo = opands2.operator()
 705     opa = opands2.operands()
 706     if opo == cot : 
 707       if filter(opa[0]) : return  -opa[0]
 708     if opo == tan : 
 709       if filter(opa[0]) : return opa[2]-pi/2
 710     return opor (*opands)
 711 
 712   if opor == arccot :
 713     opo = opands[0].operator()
 714     opa = opands[0].operands()
 715     if opo == cot : 
 716       if filter(opa[0]) : return opa[0]
 717     if opo == tan : 
 718       if filter(opa[0]) : return pi/2-opa[0] 
 719     opands2 = - opands[0]
 720     opo = opands2.operator()
 721     opa = opands2.operands()
 722     if opo == cot : 
 723       if filter(opa[0]) : return - opa[0]
 724     if opo == tan : 
 725       if filter(opa[0]) : return opa[0] - pi/2 
 726     opands2 = 1/opands[0]
 727     opo = opands2.operator()
 728     opa = opands2.operands()
 729     if opo == tan : 
 730       if filter(opa[0]) : return  opa[0]
 731     if opo == cot : 
 732       if filter(opa[0]) : return pi/2-opa[0] 
 733     opands2 = -opands2
 734     opo = opands2.operator()
 735     opa = opands2.operands()
 736     if opo == tan : 
 737       if filter(opa[0]) : return  -opa[0]
 738     if opo == cot : 
 739       if filter(opa[0]) : return opa[2]-pi/2
 740     return opor (*opands)
 741 
 742   return opor (*opands)