attachment:rewrite-20101019.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]== "asin2acos" :
 205       expr = mapexpression (expr, asin2acos, filter, level)
 206     elif rule[k]== "acos2asin" :
 207       expr = mapexpression (expr, acos2asin, filter, level)
 208     elif rule[k]== "atrigo2log" :
 209       expr = mapexpression (expr, atrigo2log, filter, level)
 210     elif rule[k]== "atrigh2log" :
 211       expr = mapexpression (expr, atrigh2log, filter, level)
 212     elif rule[k]== "factorial2gamma" :
 213       expr = mapexpression (expr, factorial2gamma, filter, level)
 214     elif rule[k]== "gamma2factorial" :
 215       expr = mapexpression (expr, gamma2factorial, filter, level)
 216 #    elif rule[k]== "binomial2factorial" :
 217 #      expr = mapexpression (expr, binomial2factorial, filter, level)
 218 
 219     elif rule[k]== "pseudoParity" :
 220       expr = mapexpression (expr, pseudoParity, filter, level)
 221 
 222     elif rule[k]== "logMainBranch" :
 223       expr = mapexpression (expr, logMainBranch, filter, level)
 224     elif rule[k]== "atrigoMainBranch" :
 225       expr = mapexpression (expr, atrigoMainBranch, filter, level)
 226 
 227 
 228     else : return "unknown rule"
 229     k=k+1
 230   return expr
 231 
 232 def positiveSin (expr) : 
 233   if pseudoPositive (expr) : return sin (expr)
 234   else : return -sin(-expr)
 235 def positiveCos (expr) : 
 236   if pseudoPositive (expr) : return cos (expr)
 237   else : return cos(-expr)
 238 
 239 def positiveTan (expr) : 
 240   if pseudoPositive (expr) : return tan (expr)
 241   else : return -tan(-expr)
 242 def positiveCot (expr) : 
 243   if pseudoPositive (expr) : return cot (expr)
 244   else : return -cot(-expr)
 245 
 246 def positiveCsc (expr) : 
 247   if pseudoPositive (expr) : return csc (expr)
 248   else : return -csc(-expr)
 249 def positiveSec (expr) : 
 250   if pseudoPositive (expr) : return sec (expr)
 251   else : return sec(-expr)
 252 
 253 def positiveSinh (expr) : 
 254   if pseudoPositive (expr) : return sinh (expr)
 255   else : return -sinh(-expr)
 256 def positiveCosh (expr) : 
 257   if pseudoPositive (expr) : return cosh (expr)
 258   else : return cosh(-expr)
 259 
 260 def positiveTanh (expr) : 
 261   if pseudoPositive (expr) : return tanh (expr)
 262   else : return -tanh(-expr)
 263 def positiveCoth (expr) : 
 264   if pseudoPositive (expr) : return coth (expr)
 265   else : return -coth(-expr)
 266 
 267 def positiveCsch (expr) : 
 268   if pseudoPositive (expr) : return csch (expr)
 269   else : return -csch(-expr)
 270 def positiveSech (expr) : 
 271   if pseudoPositive (expr) : return sech (expr)
 272   else : return sech(-expr)
 273 
 274 def positiveAsin (expr) : 
 275   if pseudoPositive (expr) : return asin (expr)
 276   else : return -asin(-expr)
 277 def positiveAtan (expr) : 
 278   if pseudoPositive (expr) : return atan (expr)
 279   else : return -atan(-expr)
 280 def positiveAcot (expr) : 
 281   if pseudoPositive (expr) : return acot (expr)
 282   else : return -acot(-expr)
 283 def positiveAcsc (expr) : 
 284   if pseudoPositive (expr) : return acsc (expr)
 285   else : return -acsc(-expr)
 286 def positiveAsinh (expr) : 
 287   if pseudoPositive (expr) : return asinh (expr)
 288   else : return -asinh(-expr)
 289 def positiveAtanh (expr) : 
 290   if pseudoPositive (expr) : return atanh (expr)
 291   else : return -atanh(-expr)
 292 def positiveAcoth (expr) : 
 293   if pseudoPositive (expr) : return acoth (expr)
 294   else : return -acoth(-expr)
 295 def positiveAcsch (expr) : 
 296   if pseudoPositive (expr) : return acsch (expr)
 297   else : return -acsch(-expr)
 298 
 299 def exp2sinhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 300   if opor == exp :
 301     if not filter(opands[0]) : return opor (*opands)
 302     if pseudoPositive (opands[0]) :
 303       return (cosh (opands[0]) + sinh (opands[0]))
 304     else : 
 305        return (cosh (-opands[0]) - sinh (-opands[0]))
 306   return opor (*opands)
 307 
 308 def exp2sincos (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 309   if opor == exp :
 310     if not filter(opands[0]) : return opor (*opands)
 311     co = I*opands[0] 
 312     if pseudoPositive (co) : return cos(co) - I * sin(co)
 313     else : return cos(-co) + I * sin(-co)
 314   return opor (*opands)
 315 
 316 def sinhcosh2exp (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 317   if opor == sinh :
 318     if not filter(opands[0]) : return opor (*opands)
 319     return (exp (opands[0]) - exp (-opands[0]))/2
 320   elif opor == cosh :
 321     if not filter(opands[0]) : return opor (*opands)
 322     return (exp (opands[0]) + exp (-opands[0]))/2
 323   return opor (*opands)
 324 
 325 def sincos2exp (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 326   if opor == sin :
 327     if not filter(opands[0]) : return opor (*opands)
 328     return -I/2 * (exp (I*opands[0]) - exp (-I*opands[0]))
 329   elif opor == cos :
 330     if not filter(opands[0]) : return opor (*opands)
 331     return (exp (I*opands[0]) + exp (-I*opands[0]))/2
 332   return opor (*opands)
 333 
 334 def trigo2sincos (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 335   if opor == tan :
 336     if not filter(opands[0]) : return opor (*opands)
 337     return sin (opands[0]) / cos (opands[0])
 338   elif opor == cot :
 339     if not filter(opands[0]) : return opor (*opands)
 340     return cos (opands[0]) / sin (opands[0])
 341   elif opor == csc :
 342     if not filter(opands[0]) : return opor (*opands)
 343     return 1 / sin (opands[0])
 344   elif opor == sec : 
 345     if not filter(opands[0]) : return opor (*opands)
 346     return 1 / cos (opands[0])
 347   return opor (*opands)
 348 
 349 def trigh2sinhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 350   if opor == tanh :
 351     if not filter(opands[0]) : return opor (*opands)
 352     return sinh (opands[0]) / cosh (opands[0])
 353   elif opor == coth :
 354     if not filter(opands[0]) : return opor (*opands)
 355     return cosh (opands[0]) / sinh (opands[0])
 356   elif opor == csch :
 357     if not filter(opands[0]) : return opor (*opands)
 358     return 1 / sinh (opands[0])
 359   elif opor == sech :
 360     if not filter(opands[0]) : return opor (*opands)
 361     return 1 / cosh (opands[0])
 362   return opor (*opands)
 363 
 364 def squareInPow(opor, opands, filter) :
 365   if (opor != operator.pow) : return(opor (*opands))
 366   if not filter[1](opands[0]) : return opor (*opands)
 367   expo = opands[1]
 368   if not (expo._is_real()) : return(opor (*opands))
 369   opo = opands[0].operator()
 370   if not opo : return(opor (*opands))
 371   if opo != filter[0] : return(opor (*opands))
 372   opa = opands[0].operands()[0]
 373   expo1 = floor(expo/2)
 374   expo2 = expo - 2*expo1
 375   return filter[2](opa)^(expo2) * filter[3](opa)^(expo1)
 376 
 377 def lessIinExp (opor, opands, filter) :
 378   if opor == exp :
 379     if not filter(opands[0]) : return opor (*opands)
 380     res = pseudoRealAndImag (opands[0]) 
 381     return exp(res[0]) * (cos(res[1]) + I*sin(res[1]))
 382   return opor (*opands)
 383 
 384 def exp2trig (opor, opands, filter) :
 385   if opor == exp :
 386     if not filter(opands[0]) : return opor (*opands)
 387     res = pseudoRealAndImag (opands[0])
 388     if pseudoPositive (res[0]) : pr = cosh(res[0])+sinh(res[0])
 389     else : pr = cosh(-res[0])-sinh(-res[0])
 390     if pseudoPositive (res[1]) : pi = cos(res[1])+I*sin(res[1]) 
 391     else : pi = cos(-res[1])-I*sin(-res[1])
 392     return  pr * pi
 393   return opor (*opands)
 394 
 395 def lessIinExp (opor, opands, filter) :
 396   if (opor == exp) and (opands[0].has(i))  :
 397     if not filter(opands[0]) : return opor (*opands)
 398     res = pseudoRealAndImag (opands[0])
 399     if pseudoPositive (res[1]) :
 400       return exp(res[0]) * (cos(res[1]) + I*sin(res[1]))
 401     else :
 402       return exp(res[0]) * (cos(-res[1]) - I*sin(-res[1]))
 403   return (opor(*opands))
 404 
 405 def lessIinTrig (opor, opands, filter) :
 406   t1 = opands[0].has(i)
 407   if not t1 :   return (opor(*opands))
 408   if opor in [sin,cos,tan,cot,csc,sec,sinh,cosh,tanh,coth,csch,sech] :
 409     if not filter(opands[0]) : return opor (*opands)
 410     res = pseudoRealAndImag (opands[0]) 
 411     if res[0] != 0 : return (opor(*opands))
 412   if opor == sin : return i*positiveSin(res[1])
 413   elif opor == cos : return positiveCosh(res[1])
 414   elif opor == tan : return i*positiveTanh(res[1])
 415   elif opor == cot : return -i*positiveCoth(res[1])
 416   elif opor == sec : return positiveSech(res[1])
 417   elif opor == csc : return -i/positiveCsch(res[1])
 418   elif opor == sinh : return i*positiveSin(res[1])
 419   elif opor == cosh : return positiveCos(res[1])
 420   elif opor == tanh : return i*positiveTan(res[1])
 421   elif opor == coth : return -i*positiveCot(res[1])
 422   elif opor == sech : return positiveSech(res[1])
 423   elif opor == csch : return -i/positiveCsc(res[1])
 424   return (opor(*opands))
 425 
 426 def sin2tancos (opor, opands, filter)  : 
 427   if opor == sin :
 428     if not filter(opands[0]) : return opor (*opands)
 429     co = opands[0] 
 430     if pseudoPositive (co) : return cos(co) * tan(co)
 431     else : return - cos(-co) * tan(-co)
 432   return opor (*opands)
 433 
 434 def sinh2tanhcosh (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 435   if opor == sinh :
 436     if not filter(opands[0]) : return opor (*opands)
 437     co = opands[0] 
 438     if pseudoPositive (co) : return cosh(co) * tanh(co)
 439     else : return - cosh(-co) * tanh(-co)
 440   return opor (*opands)
 441 
 442 def sincos2tanHalf (opor, opands, filter)  : 
 443   co = opands[0]/2 
 444   if opor == sin :
 445     if not filter(opands[0]) : return opor (*opands)
 446     if pseudoPositive (co/2) : return 2*tan(co)/(1+tan(co)^2)
 447     else : return 2*tan(-co)/(1+tan(-c)^2)
 448   if opor == cos :
 449     if pseudoPositive (co/2) : return (1-tan(co)^2)/(1+tan(co)^2)
 450     else : return (1-tan(-co)^2)/(1+tan(-c)^2)
 451   return opor (*opands)
 452 
 453 def sinhcosh2tanhHalf (opor, opands, filter)  : # exp(x) => sinh(x) + cosh(x)
 454   co = opands[0]/2 
 455   if opor == sinh :
 456     if not filter(opands[0]) : return opor (*opands)
 457     if pseudoPositive (co/2) : return 2*tanh(co)/(1-tanh(co)^2)
 458     else : return 2*tanh(-co)/(1-tanh(-c)^2)
 459   if opor == cosh :
 460     if not filter(opands[0]) : return opor (*opands)
 461     if pseudoPositive (co/2) : return (1+tanh(co)^2)/(1-tanh(co)^2)
 462     else : return (1+tanh(-co)^2)/(1-tanh(-c)^2)
 463   return opor (*opands)
 464 
 465 def pseudoParity (opor, opands, filter) :
 466   if opor == sin : 
 467     if not filter(opands[0]) : return opor (*opands)
 468     return positiveSin (opands[0])
 469   elif opor == cos :
 470     if not filter(opands[0]) : return opor (*opands)
 471     return positiveCos (opands[0])
 472   elif opor == tan :
 473     if not filter(opands[0]) : return opor (*opands)
 474     return positiveTan (opands[0])  
 475   elif opor == cot :
 476     if not filter(opands[0]) : return opor (*opands)
 477     return positiveCot (opands[0])  
 478   elif opor == sec :
 479     if not filter(opands[0]) : return opor (*opands)
 480     return positiveSec (opands[0])  
 481   elif opor == csc :
 482     if not filter(opands[0]) : return opor (*opands)
 483     return positiveCsc (opands[0])
 484 
 485   elif opor == sinh :
 486     if not filter(opands[0]) : return opor (*opands)
 487     return positiveSinh (opands[0])
 488   elif opor == cosh :
 489     if not filter(opands[0]) : return opor (*opands)
 490     return positiveCosh (opands[0])
 491   elif opor == tanh :
 492     if not filter(opands[0]) : return opor (*opands)
 493     return positiveTanh (opands[0])  
 494   elif opor == coth :
 495     if not filter(opands[0]) : return opor (*opands)
 496     return positiveCoth (opands[0])  
 497   elif opor == sech :
 498     if not filter(opands[0]) : return opor (*opands)
 499     return positiveSech (opands[0])  
 500   elif opor == csch :
 501     if not filter(opands[0]) : return opor (*opands)
 502     return positiveCsch (opands[0])
 503 
 504   elif opor == asin :
 505     if not filter(opands[0]) : return opor (*opands)
 506     return positiveAsin (opands[0])
 507   elif opor == atan :
 508     if not filter(opands[0]) : return opor (*opands)
 509     return positiveAtan (opands[0])  
 510   elif opor == acot :
 511     if not filter(opands[0]) : return opor (*opands)
 512     return positiveAcot (opands[0])  
 513   elif opor == acsc :
 514     if not filter(opands[0]) : return opor (*opands)
 515     return positiveAcsc (opands[0])  
 516 
 517   elif opor == asinh :
 518     if not filter(opands[0]) : return opor (*opands)
 519     return positiveAsinh (opands[0])
 520   elif opor == atanh :
 521     if not filter(opands[0]) : return opor (*opands)
 522     return positiveAtanh (opands[0])  
 523   elif opor == acoth :
 524     if not filter(opands[0]) : return opor (*opands)
 525     return positiveAcoth (opands[0])  
 526   elif opor == acsch :
 527     if not filter(opands[0]) : return opor (*opands)
 528     return positiveAcsch (opands[0])  
 529 
 530   else : return opor(*opands)
 531 
 532 def pow2exp (opor, opands, filter) :
 533   if opor == operator.pow : 
 534     if not filter(opands[0], opands[1]) : return opor (*opands)
 535     return (exp (log(opands[0]) * opands[1]))
 536   else : return opor(*opands)
 537 
 538 def exp2pow (opor, opands, filter) :
 539   if opor == exp : 
 540     if not filter(opands[0]) : return opor (*opands)
 541     res = opands[0].match(SR.wild(0)*log(SR.wild(1)))
 542     if res : return res[SR.wild(1)]^res[SR.wild(0)]
 543   return opor(*opands)
 544 
 545 def asin2acos (opor, opands, filter) :
 546   if opor == asin : 
 547     if not filter(opands[0]) : return opor (*opands)
 548     return pi/2-acos(opands[0])
 549   return opor(*opands)
 550 
 551 def acos2asin (opor, opands, filter) :
 552   if opor == acos : 
 553     if not filter(opands[0]) : return opor (*opands)
 554     return pi/2-asin(opands[0])
 555   return opor(*opands)
 556 
 557 def atrigh2log (opor, opands, filter) :
 558   if opor == asinh : 
 559     if not filter(opands[0]) : return opor (*opands)
 560     return log(opands[0]+sqrt(opands[0]^2+1))
 561   elif opor == acosh : 
 562     if not filter(opands[0]) : return opor (*opands)
 563     return log(opands[0]+sqrt(opands[0]^2-1))
 564   elif opor == atanh :
 565     if not filter(opands[0]) : return opor (*opands)
 566     return 1/2*(log(1+opands[0])-log(1-opands[0]))
 567   elif opor == acoth :
 568     if not filter(opands[0]) : return opor (*opands)
 569     return 1/2*(log(1+1/opands[0])-log(1-1/opands[0]))
 570   elif opor == acsch : 
 571     if not filter(opands[0]) : return opor (*opands)
 572     return log(1/opands[0]+sqrt(1/opands[0]^2+1)) 
 573   elif opor == asech : 
 574     if not filter(opands[0]) : return opor (*opands)
 575     return log(1/opands[0]+sqrt(1/opands[0]^2-1))
 576   return opor(*opands)
 577 
 578 def atrigo2log (opor, opands, filter) :
 579   if opor == asin : 
 580     if not filter(opands[0]) : return opor (*opands)
 581     return -i*log(i*opands[0]+sqrt(1-opands[0]^2))
 582   elif opor == acos : 
 583     if not filter(opands[0]) : return opor (*opands)
 584     return -i*log(opands[0]+i*sqrt(1-opands[0]^2))
 585   elif opor == atan : 
 586     if not filter(opands[0]) : return opor (*opands)
 587     return i/2*log((1-i*opands[0])-log(1+i*opands[0]))
 588   elif opor == acot : 
 589     if not filter(opands[0]) : return opor (*opands)
 590     return i/2*log((1-i/opands[0])-(1+i/opands[0]))
 591   elif opor == acsc : 
 592     if not filter(opands[0]) : return opor (*opands)
 593     return -i*log(i/opands[0]+sqrt(1-1/opands[0]^2))
 594   elif opor == asec : 
 595     if not filter(opands[0]) : return opor (*opands)
 596     return -i*log(1/opands[0]+i*sqrt(1-1/opands[0]^2))
 597 
 598 def factorial2gamma (opor, opands, filter) :
 599   if opor == factorial : 
 600     if not filter(opands[0]) : return opor (*opands)
 601     return gamma(opands[0]+1)
 602   return opor(*opands)
 603 def gamma2factorial (opor, opands, filter) :
 604   if opor == gamma : 
 605     if not filter(opands[0]) : return opor (*opands)
 606     return factorial(opands[0]-1)
 607   return opor(*opands)
 608 
 609 #def binomial2factorial (opor, opands, filter) :
 610 #  if opor == binomial : 
 611 #    return    factorial (opand[0]) \
 612 #           / (factorial (operands[1]) * factorial (operands[0]-operands[1]))
 613 #  return opor(*opands)
 614 
 615 def logMainBranch (opor, opands, filter) :
 616 #  if opor != log : return(opor (*opands))
 617   if opor != ln : return(opor (*opands))
 618   opo = opands[0].operator()
 619   opa = opands[0].operands()
 620   if opo != exp : return(opor (*opands))
 621   if not filter(opa[0]) : return opor (*opands)
 622   else : return opa[0]
 623   return(opor (*opands))
 624 
 625 # atrigBranchCut : [acos|asin] o [+|-] [sin|cos] = 8 cas
 626 #                  [atan|acot] o [+|-] o [id|inverse] [tan|cot] = 16 cas
 627 
 628 def atrigoMainBranch (opor, opands, filter) :
 629   if opor == acos :
 630     opo = opands[0].operator()
 631     opa = opands[0].operands()
 632     if opo == cos : 
 633       if filter(opa[0]) : return opa[0]
 634     if opo == sin : 
 635       if filter(opa[0]) : return pi/2-opa[0] 
 636     opands2 = - opands[0]
 637     opo = opands2.operator()
 638     opa = opands2.operands()
 639     if opo == cos : 
 640       if filter(opa[0]) : return pi - opa[0]
 641     if opo == sin : 
 642       if filter(opa[0]) : return pi/2 - opa[0]
 643     return opor (*opands)
 644 
 645   if opor == asin :
 646     opo = opands[0].operator()
 647     opa = opands[0].operands()
 648     if opo == sin : 
 649       if filter(opa[0]) : return opa[0]
 650     if opo == cos : 
 651       if filter(opa[0]) : return pi/2 - opa[0] 
 652     opands2 = - opands[0]
 653     opo = opands2.operator()
 654     opa = opands2.operands()
 655     if opo == sin : 
 656       if filter(opa[0]) : return - opa[0]
 657     if opo == cos : 
 658       if filter(opa[0]) : return opa[0] - pi/2
 659     return opor (*opands)
 660 
 661   if opor == atan :
 662     opo = opands[0].operator()
 663     opa = opands[0].operands()
 664     if opo == tan : 
 665       if filter(opa[0]) : return opa[0]
 666     if opo == cot : 
 667       if filter(opa[0]) : return pi/2-opa[0] 
 668     opands2 = - opands[0]
 669     opo = opands2.operator()
 670     opa = opands2.operands()
 671     if opo == tan : 
 672       if filter(opa[0]) : return - opa[0]
 673     if opo == cot : 
 674       if filter(opa[0]) : return opa[0] - pi/2 
 675     opands2 = 1/opands[0]
 676     opo = opands2.operator()
 677     opa = opands2.operands()
 678     if opo == cot : 
 679       if filter(opa[0]) : return  opa[0]
 680     if opo == tan : 
 681       if filter(opa[0]) : return pi/2-opa[0] 
 682     opands2 = -opands2
 683     opo = opands2.operator()
 684     opa = opands2.operands()
 685     if opo == cot : 
 686       if filter(opa[0]) : return  -opa[0]
 687     if opo == tan : 
 688       if filter(opa[0]) : return opa[2]-pi/2
 689     return opor (*opands)
 690 
 691   if opor == acot :
 692     opo = opands[0].operator()
 693     opa = opands[0].operands()
 694     if opo == cot : 
 695       if filter(opa[0]) : return opa[0]
 696     if opo == tan : 
 697       if filter(opa[0]) : return pi/2-opa[0] 
 698     opands2 = - opands[0]
 699     opo = opands2.operator()
 700     opa = opands2.operands()
 701     if opo == cot : 
 702       if filter(opa[0]) : return - opa[0]
 703     if opo == tan : 
 704       if filter(opa[0]) : return opa[0] - pi/2 
 705     opands2 = 1/opands[0]
 706     opo = opands2.operator()
 707     opa = opands2.operands()
 708     if opo == tan : 
 709       if filter(opa[0]) : return  opa[0]
 710     if opo == cot : 
 711       if filter(opa[0]) : return pi/2-opa[0] 
 712     opands2 = -opands2
 713     opo = opands2.operator()
 714     opa = opands2.operands()
 715     if opo == tan : 
 716       if filter(opa[0]) : return  -opa[0]
 717     if opo == cot : 
 718       if filter(opa[0]) : return opa[2]-pi/2
 719     return opor (*opands)
 720 
 721   return opor (*opands)