def test_asinh(self):
        import math

        self.ftest(math.asinh(0), 0)
        self.ftest(math.asinh(1), 0.88137358701954305)
        self.ftest(math.asinh(-1), -0.88137358701954305)
        assert math.isinf(math.asinh(float("inf")))

 def testAsinh(self):
     self.assertRaises(TypeError, math.asinh)
     self.ftest('asinh(0)', math.asinh(0), 0)
     self.ftest('asinh(1)', math.asinh(1), 0.88137358701954305)
     self.ftest('asinh(-1)', math.asinh(-1), -0.88137358701954305)
     self.assertEqual(math.asinh(INF), INF)
     # self.assertEqual(math.asinh(NINF), NINF)
     self.assertTrue(math.isnan(math.asinh(NAN)))

def g(p):
  x, y, z = p[0], p[1], abs(p[2])
  return + x*y*z * asinh(z / (sqrt(x**2 + y**2) + eps))                         \
         + y / 6.0 * (3.0 * z**2 - y**2) * asinh(x / (sqrt(y**2 + z**2) + eps)) \
         + x / 6.0 * (3.0 * z**2 - x**2) * asinh(y / (sqrt(x**2 + z**2) + eps)) \
         - z**3 / 6.0 * atan(x*y / (z * sqrt(x**2 + y**2 + z**2) + eps))        \
         - z * y**2 / 2.0 * atan(x*z / (y * sqrt(x**2 + y**2 + z**2) + eps))    \
         - z * x**2 / 2.0 * atan(y*z / (x * sqrt(x**2 + y**2 + z**2) + eps))    \
         - x*y * sqrt(x**2 + y**2 + z**2) / 3.0

def tri_angle_side_angle(a1, s, a2):
  """if the sides go [s,X,Y] with angles a1, a2, on either side of s, this 
  function returns X,Y"""
  #find the remaining angle
  S = math.acos(math.sin(a1)*math.sin(a2)*math.cosh(s) - math.cos(a1)*math.cos(a2))
  #use the law of sines to get the other sides
  sin_ratio = math.sinh(s)/math.sin(S)
  X = math.asinh(sin_ratio*math.sin(a1))
  Y = math.asinh(sin_ratio*math.sin(a2))
  return X,Y

def coefRAsientoZapRectangularFlexible(l,b):
    '''
    coefRAsientoZapRectangularFlexible(l,b):
    Devuelve el coeficiente R de la expresión del asiento
    de una cimentación superficial flexible rectangular de
    acuerdo con el apartado 4.8.1 de la Guia de
    cimentaciones en obras de carretera, figura 4.10 de la página 113
    '''
    retv=b*math.asinh(l/float(b))+l*math.asinh(b/float(l))
    return(retv)

def acos(x):
    _acos_special = [
        [3*pi/4+infj, pi+infj, pi+infj, pi-infj, pi-infj, 3*pi/4-infj, nan+infj],
        [pi/2+infj, None, None, None, None, pi/2-infj, nan+nanj],
        [pi/2+infj, None, None, None, None, pi/2-infj, pi/2+nanj],
        [pi/2+infj, None, None, None, None, pi/2-infj, pi/2+nanj],
        [pi/2+infj, None, None, None, None, pi/2-infj, nan+nanj],
        [pi/4+infj, infj, infj, 0.0-infj, 0.0-infj, pi/4-infj, nan+infj],
        [nan+infj, nan+nanj, nan+nanj, nan+nanj, nan+nanj, nan-infj, nan+nanj]
    ]

    z = _make_complex(x)

    if not isfinite(z):
        return _acos_special[_special_type(z.real)][_special_type(z.imag)]

    if abs(z.real) > _LARGE_DOUBLE or abs(z.imag) > _LARGE_DOUBLE:
        if z.real < 0:
            imag = -math.copysign(math.log(math.hypot(z.real/2, z.imag/2)) +
                                  2 * _LOG_2, z.imag)
        else:
            imag = math.copysign(math.log(math.hypot(z.real/2, z.imag/2)) +
                                 2 * _LOG_2, -z.imag)
        return complex(math.atan2(abs(z.imag), z.real), imag)

    s1 = sqrt(complex(1.0 - z.real, -z.imag))
    s2 = sqrt(complex(1.0 + z.real, z.imag))
    return complex(2 * math.atan2(s1.real, s2.real),
                   math.asinh(s2.real*s1.imag - s2.imag*s1.real))

    def test_acosh(self):
        import math

        self.ftest(math.acosh(1), 0)
        self.ftest(math.acosh(2), 1.3169578969248168)
        assert math.isinf(math.asinh(float("inf")))
        raises(ValueError, math.acosh, 0)

def deg2num(x, y, zoom):
    """Convert longitude, latitude to tile numbers."""
    xmerc = math.radians(x)
    ymerc = math.asinh(math.tan(math.radians(y)))
    xtile = int((1 + xmerc/math.pi) / 2 * 2**zoom)
    ytile = int((1 - ymerc/math.pi) / 2 * 2**zoom)
    return xtile, ytile

  def zetainv(self, taup, lam):
    psi = math.asinh(taup)
    
    scal = 1.0 / math.hypot(1.0, taup);

    (u, v, flag) = self.zetainv0(psi, lam)
    if flag:
      return [u, v]

    stol2 = self.tolerance2 / max(psi, 1.0)**2.0
    # min iterations = 2, max iterations = 6; mean = 4.0
    trip = 0
    
    for i in range(100):
      (snu, cnu, dnu, ph) = ellipj(u, self.mu)
      (snv, cnv, dnv, ph) = ellipj(v, self.mv)
      (tau1, lam1) =  self.zeta(snu, cnu, dnu, snv, cnv, dnv);
      (du1, dv1) = self.dwdzeta(snu, cnu, dnu, snv, cnv, dnv)
      
      tau1 -= taup;  lam1 -= lam;
      tau1 *= scal;
      delu = tau1 * du1 - lam1 * dv1
      delv = tau1 * dv1 + lam1 * du1
      u -= delu;  v -= delv;
      if trip > 0:
        break;
      delw2 = delu*delu + delv*delv;
      if delw2 < stol2 :
        trip += 1
    return [u, v]

 def zeta(self, snu, cnu, dnu, snv, cnv, dnv):
   '''
    Lee 54.17 but write
    atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2))
    atanh(_e * snu / dnv) =
            asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2))
   '''
   d1 = math.sqrt(cnu*cnu + self.mv * (snu*snv)**2.0)
   d2 = math.sqrt(self.mu * cnu*cnu + self.mv * cnv*cnv)
   t1 = snu * dnv
   if d1 != 0.0:
     t1 /= d1
   else: 
     t1 = self.overflow if snu >= 0.0 else -self.overflow
   
   t2 = self.e * snu
   if d2 != 0.0:
     t2 = math.sinh(self.e * math.asinh(t2 / d2))
   else:
     t2 = self.overflow if snu >= 0.0 else -self.overflow
   
   taup = t1 * math.hypot(1.0, t2) - t2 * math.hypot(1.0, t1);
   lam = 0.0
   if d1 != 0.0 and d2 != 0.0:
     lam = math.atan2(dnu*snv, cnu*cnv) - self.e * math.atan2(self.e*cnu*snv, dnu*cnv)
   return [taup, lam]

def hyperbolic_fun_cal(inp_val1, opn_type):
     oprn_dic = {
        '1': 'inverse hyperbolic cosine of x', '2':'inverse hyperbolic sine of x',
        '3':' inverse hyperbolic tangent of x', '4':'hyperbolic cosine of x',
        '5':'hyperbolic sine of x', '6':'hyperbolic tangent of x'}
     if int(opn_type) == 1:
         output = math.acosh(float(inp_val1))
         return str(output)
     if int(opn_type) == 2:
         output = math.asinh(float(inp_val1))
         return str(output)
     if int(opn_type) == 3:
         output = math.atanh(float(inp_val1))
         return str(output)
     if int(opn_type) == 4:
         output = math.cosh(float(inp_val1))
         return str(output)
     if int(opn_type) == 5:
         output = math.sinh(float(inp_val1))
         return str(output)
     if int(opn_type) == 6:
         output = math.tanh(float(inp_val1))
         return str(output)
     else:
         return "Invalid Operation"

def acos(x):
    """
        Return the arc cosine of x.

        There are two branch cuts: One extends right from 1 along the real axis to ∞, continuous from below.
        The other extends left from -1 along the real axis to -∞, continuous from above.
    """

    ret = _SPECIAL_VALUE(x, _acos_special_values)
    if ret is not None:
        return ret

    if math.fabs(x.real) > _CM_LARGE_DOUBLE or math.fabs(x.imag) > _CM_LARGE_DOUBLE:

        # avoid unnecessary overflow for large arguments
        _real = math.atan2(math.fabs(x.imag), x.real)

        # split into cases to make sure that the branch cut has the
        # correct continuity on systems with unsigned zeros
        if x.real < 0:
            _imag = -math.copysign(math.log(math.hypot(x.real/2., x.imag/2.)) + _M_LN2*2., x.imag);
        else:
            _imag = math.copysign(math.log(math.hypot(x.real/2., x.imag/2.)) + _M_LN2*2., -x.imag);
    else:
        s1 = complex(float(1-x.real), -x.imag)
        s1 = sqrt(s1)
        s2 = complex(1.0+x.real, x.imag)
        s2 = sqrt(s2)
        _real = 2.0*math.atan2(s1.real, s2.real);
        _imag = math.asinh(s2.real*s1.imag - s2.imag*s1.real)

    return complex(_real,_imag)

def asinh(x):
    _asinh_special = [
        [-inf-1j*pi/4, complex(-float("inf"), -0.0), complex(-float("inf"), -0.0),
            complex(-float("inf"), 0.0), complex(-float("inf"), 0.0), -inf+1j*pi/4, -inf+nanj],
        [-inf-1j*pi/2, None, None, None, None, -inf+1j*pi/2, nan+nanj],
        [-inf-1j*pi/2, None, None, None, None, -inf+1j*pi/2, nan+nanj],
        [inf-1j*pi/2, None, None, None, None, inf+1j*pi/2, nan+nanj],
        [inf-1j*pi/2, None, None, None, None, inf+1j*pi/2, nan+nanj],
        [inf-1j*pi/4, complex(float("inf"), -0.0), complex(float("inf"), -0.0),
            inf, inf, inf+1j*pi/4, inf+nanj],
        [inf+nanj, nan+nanj, complex(float("nan"), -0.0), nan, nan+nanj, inf+nanj, nan+nanj]
    ]

    z = _make_complex(x)

    if not isfinite(z):
        return _asinh_special[_special_type(z.real)][_special_type(z.imag)]

    if abs(z.real) > _LARGE_DOUBLE or abs(z.imag) > _LARGE_DOUBLE:
        if z.imag >= 0:
            real = math.copysign(math.log(math.hypot(z.imag/2, z.real/2)) +
                                 2 * _LOG_2, z.real)
        else:
            real = -math.copysign(math.log(math.hypot(z.imag/2, z.real/2)) +
                                  2 * _LOG_2, -z.real)
        return complex(real, math.atan2(z.imag, abs(z.real)))

    s1 = sqrt(complex(1+z.imag, -z.real))
    s2 = sqrt(complex(1-z.imag, z.real))
    return complex(math.asinh(s1.real*s2.imag-s2.real*s1.imag),
                   math.atan2(z.imag, s1.real*s2.real - s1.imag*s2.imag))

    def eccentric_anomaly(self, time):
        """Eccentric anomaly at a given time (s)"""
        M = self.mean_anomaly(time)
        e = self.eccentricity

        if e < 1:  # M = E - e sin E
            M %= (2*math.pi)

            # sin(E) = E -> M = (1 - e) E
            if abs(M) < 2**-26:
                return M / (1 - e)

            return analysis.newton_raphson(
                x_0=math.pi,
                f=lambda E: E - e*math.sin(E) - M,
                f_prime=lambda E: 1 - e*math.cos(E),
            )
        else:  # M = e sinh E - E
            # sinh(E) = E -> M = (e - 1) E
            if abs(M) < 2**-26:
                return M / (e - 1)

            return analysis.newton_raphson(
                x_0=math.asinh(M),
                f=lambda E: e*math.sinh(E) - E - M,
                f_prime=lambda E: e*math.cosh(E) - 1,
            )

    def sinh_effective_mass_errors(self, dt, fast=True, period=None):
        if fast: logging.info("sinh emass computed fast method")

        T = self.period_check(period)

        jkasv = self.jackknife_average_sub_vev()
        jkemass = {}
        for cfg in self.configs:
            asvc = jkasv[cfg]
            emass = {}
            for t in self.times[dt:-dt]:
                if t in self.emass_skip_times:
                    emass[t] = 0.0
                    continue
                try:
                    guess = (1.0 / float(dt))*math.asinh((asvc[t+dt] + asvc[t-dt])/(2.0*asvc[t]))
                    if fast:
                        emass[t] = guess
                    else:
                        emass[t] = newton.newton_sinh_for_m(t,t+dt,asvc, guess, T)
                except ValueError:
                    #logging.debug("invalid argument to log, setting to zero")
                    emass[t] = 0.0
                except ZeroDivisionError:
                    logging.debug("div by zero, setting to zero")
                    emass[t] = 0.0
                except KeyError:
                    logging.error("index out of range")
            jkemass[cfg] = emass
        jkemassobj = configtimeobj.Cfgtimeobj.fromDataDict(jkemass)
        effmass_dt = self.sinh_effective_mass(dt, fast=fast, period=T)
        return {t: jackknife.errorbars(effmass_dt[t], jkemassobj.get(time=t))
                for t in self.times[dt:-dt]}

    def test_arcsinh(self):
        import math
        from _numpypy import arcsinh

        for v in [float('inf'), float('-inf'), 1.0, math.e]:
            assert math.asinh(v) == arcsinh(v)
        assert math.isnan(arcsinh(float("nan")))

def pmuboost3d(jets, met, lep):
    mw = 80.385
    wwwtlv = ROOT.TLorentzVector(1e-9,1e-9,1e-9,1e-9)
    mettlv = ROOT.TLorentzVector(1e-9,1e-9,1e-9,1e-9)
    leptlv = ROOT.TLorentzVector(1e-9,1e-9,1e-9,1e-9)
    for jet in jets:
#        if(e.jetPt[ij] > 30.):
            tlvaux = ROOT.TLorentzVector(1e-9,1e-9,1e-9,1e-9)
            tlvaux.SetPtEtaPhiM(jet['pt'],0.,jet['phi'],0.)
            wwwtlv -= tlvaux
    mettlv.SetPtEtaPhiM(met['pt'],0.,met['phi'],0.)
    leptlv.SetPtEtaPhiM(lep['pt'],lep['eta'],lep['phi'], 0.)
#    leptlv.SetPtEtaPhiM(e.muPt[0],e.muEta[0],e.muPhi[0],0.)
# calculate eta(W) estimate; take the + solution
    ptl = (wwwtlv - mettlv).Pt()
    ptn = mettlv.Pt()
    ptw = wwwtlv.Pt()
    etal = leptlv.Eta()
    A = (mw*mw + ptw*ptw - ptl*ptl - ptn*ptn) / (2.*ptl*ptn)
    if(A<1.): A=1.
    etaw = asinh((ptl*sinh(etal)+ptn*sinh(etal+acosh(A)))/ptw)
    if(abs(etaw)>10.): etaw=10*etaw/abs(etaw) # avoid too large values
    wwwtlv.SetPtEtaPhiM(wwwtlv.Pt(),etaw,wwwtlv.Phi(),mw)
# calulate boosted lepton momentum
    beta = wwwtlv.BoostVector()
    leptlv.Boost(-beta)
    return leptlv.P()

  def test_funcs_multi(self):
    pi = math.pi

    # sin family
    self.assertQuantity(data.sin(Quantity( (0,pi/2), (2,2))), (0,1), (2,0)) 
    self.assertQuantity(data.sinh(Quantity((0,1), (2,2))), (0, math.sinh(1)), (2, math.cosh(1)*2)) 
    self.assertQuantity(data.asin(Quantity((0,0.5), (2,2))), (0,math.asin(0.5)), (2,2/math.sqrt(1-0.5**2))) 
    self.assertQuantity(data.asinh(Quantity((0,1), (2,2))), (0,math.asinh(1)), (2,2/math.sqrt(1+1**2))) 

    # cos family
    self.assertQuantity(data.cos(Quantity((0,pi/2), (2,2))), (1,0), (0,-2)) 
    self.assertQuantity(data.cosh(Quantity((0,1), (2,2))), (1,math.cosh(1)), (0,math.sinh(1)*2)) 
    self.assertQuantity(data.acos(Quantity((0,0.5), (2,2))), (math.acos(0),math.acos(0.5)), (-2,-2/math.sqrt(1-0.5**2)))
    self.assertQuantity(data.acosh(Quantity((2,3), (2,2))), (math.acosh(2), math.acosh(3)), (2/math.sqrt(2**2-1),2/math.sqrt(3**2-1)))

    # tan family
    self.assertQuantity(data.tan(Quantity((0,1), (2,2))), (0,math.tan(1)), (2,2/math.cos(1)**2))
    self.assertQuantity(data.tanh(Quantity((0,1), (2,2))), (0,math.tanh(1)), (2,2/math.cosh(1)**2)) 
    self.assertQuantity(data.atan(Quantity((0,1), (2,2))), (0, math.atan(1)), (2,2/(1+1**2))) 
    self.assertQuantity(data.atan2(Quantity((0,1), (2,2)), Quantity((1,1), (0,0))), (0,math.atan(1)), (2,2/(1+1**2))) 
    self.assertQuantity(data.atanh(Quantity((0,0.5), (2,2))), (0,math.atanh(0.5)), (2,2/(1-0.5**2)))

    #misc
    self.assertQuantity(data.sqrt(Quantity((1,4), (2,2))), (1,2), (1,1/2))
    self.assertQuantity(data.exp(Quantity((1,4), (2,2))), (math.e, math.e**4), (2 * math.e,2*math.e**4))
    self.assertQuantity(data.log(Quantity((1,4), (2,2))), (0, math.log(4)), (2,1/2))

    def compute(self, plug, data):

        #   Check if output value is connected
        if plug == self.aOutputVaue:

            #    Get input datas
            operationTypeHandle = data.inputValue(self.aOperationType)
            operationType = operationTypeHandle.asInt()

            inputValueHandle = data.inputValue(self.aInputValue)
            inputValue = inputValueHandle.asFloat()
            
            #   Math Sinus
            outputValue = 0
            if operationType == 0:
                outputValue = math.asin(inputValue)
            if operationType == 1:
                outputValue = math.sin(inputValue)
            if operationType == 2:
                outputValue = math.asinh(inputValue)
            if operationType == 3:
                outputValue = math.sinh(inputValue)

            #   Output Value
            output_data = data.outputValue(self.aOutputVaue)
            output_data.setFloat(outputValue)

        #   Clean plug
        data.setClean(plug)

    def test_arcsinh(self):
        import math
        from numpypy import arcsinh, inf

        for v in [inf, -inf, 1.0, math.e]:
            assert math.asinh(v) == arcsinh(v)
        assert math.isnan(arcsinh(float("nan")))