Python cmath.sqrt() Examples

def quadratic(a, b, c=None): 
    """ 
    x^2 + ax + b = 0 (or ax^2 + bx + c = 0) By substituting x = y-t and t = a/2,
    the equation reduces to y^2 + (b-t^2) = 0 which has easy solution 
    y = +/- sqrt(t^2-b)

    Author: Victor Gil
    """ 
    if c: # (ax^2 + bx + c = 0) 
        a, b = b / float(a), c / float(a)
    t = a / 2.0 
    r = t**2 - b 
    if r >= 0: # real roots 
        y1 = math.sqrt(r) 
    else: # complex roots 
        y1 = cmath.sqrt(r) 
    y2 = -y1 
    return y1 - t, y2 - t 

def _SA_partial_horiz_torispherical_head_int_1(x, b, c):
    x0 = x*x
    x1 = b - x0
    x2 = x1**0.5
    x3 = -b + x0
    x4 = c*c
    try:
        x5 = (x1 - x4)**-0.5
    except:
        x5 = (x1 - x4+0j)**-0.5
    x6 = x3 + x4
    x7 = b**0.5
    try:
        x3_pow = x3**(-1.5)
    except:
        x3_pow = (x3+0j)**(-1.5)
    ans = (x*cacos(c/x2) + x3_pow*x5*(-c*x1*csqrt(-x6*x6)*catan(x*x2/(csqrt(x3)*csqrt(x6)))
        + x6*x7*csqrt(-x1*x1)*catan(c*x*x5/x7))/csqrt(-x6/x1))
    return ans.real 

def transversity_amps_qcdf(q2, wc, par, B, V, **kwargs):
    """QCD factorization corrections to B->Vll transversity amplitudes."""
    mB = par['m_'+B]
    mV = par['m_'+V]
    scale = config['renormalization scale']['bvll']
    # using the b quark pole mass here!
    mb = running.get_mb_pole(par)
    N = flavio.physics.bdecays.bvll.amplitudes.prefactor(q2, par, B, V)/4
    T_perp_ = T_perp(q2, par, wc, B, V, scale, **kwargs)
    T_para_ = T_para(q2, par, wc, B, V, scale, **kwargs)
    ta = {}
    ta['perp_L'] = N * sqrt(2)*2 * (mB**2-q2) * mb / q2 * T_perp_
    ta['perp_R'] =  ta['perp_L']
    ta['para_L'] = -ta['perp_L']
    ta['para_R'] =  ta['para_L']
    ta['0_L'] = ( N * mb * (mB**2 - q2)**2 )/(mB**2 * mV * sqrt(q2)) * T_para_
    ta['0_R'] = ta['0_L']
    ta['t'] = 0
    ta['S'] = 0
    return ta 

def roots_cubic_a2(a, b, c, d):
    t2 = a*a
    t3 = d*d
    t10 = c*c
    t14 = b*b
    t15 = t14*b
    t20 = csqrt(-18.0*a*b*c*d + 4.0*a*t10*c + 4.0*t15*d - t14*t10 + 27.0*t2*t3)
    t31 = (36.0*c*b*a + 12.0*root_three*t20*a - 108.0*d*t2 - 8.0*t15)**third
    t32 = 1.0/a
    root1 = t31*t32*sixth - two_thirds*(3.0*a*c - t14)*t32/t31 - b*t32*third
    t33 = t31*t32
    t40 = (3.0*a*c - t14)*t32/t31
    
    t50 = -t33*twelfth + t40*third - b*t32*third
    t51 = 0.5j*root_three *(t33*sixth + two_thirds*t40)
    root2 = t50 + t51
    root3 = t50 - t51
    return [root1, root2, root3] 

def roots_cubic_a1(b, c, d):
    t1 = b*b
    t2 = t1*b
    t4 = c*b
    t9 = c*c
    t16 = d*d
    t19 = csqrt(12.0*t9*c + 12.0*t2*d - 54.0*t4*d - 3.0*t1*t9 + 81.0*t16)
    t22 = (-8.0*t2 + 36.0*t4 - 108.0*d + 12.0*t19)**third
    root1 = t22*sixth - 6.0*(c*third - t1*ninth)/t22 - b*third
    t28 = (c*third - t1*ninth)/t22
    t101 = -t22*twelfth + 3.0*t28 - b*third
    t102 =  root_three*(t22*sixth + 6.0*t28)
    
    root2 = t101 + 0.5j*t102
    root3 = t101 - 0.5j*t102
    
    return [root1, root2, root3] 

def test_cyclotomic():
    mp.dps = 15
    assert [cyclotomic(n,1) for n in range(31)] == [1,0,2,3,2,5,1,7,2,3,1,11,1,13,1,1,2,17,1,19,1,1,1,23,1,5,1,3,1,29,1]
    assert [cyclotomic(n,-1) for n in range(31)] == [1,-2,0,1,2,1,3,1,2,1,5,1,1,1,7,1,2,1,3,1,1,1,11,1,1,1,13,1,1,1,1]
    assert [cyclotomic(n,j) for n in range(21)] == [1,-1+j,1+j,j,0,1,-j,j,2,-j,1,j,3,1,-j,1,2,1,j,j,5]
    assert [cyclotomic(n,-j) for n in range(21)] == [1,-1-j,1-j,-j,0,1,j,-j,2,j,1,-j,3,1,j,1,2,1,-j,-j,5]
    assert cyclotomic(1624,j) == 1
    assert cyclotomic(33600,j) == 1
    u = sqrt(j, prec=500)
    assert cyclotomic(8, u).ae(0)
    assert cyclotomic(30, u).ae(5.8284271247461900976)
    assert cyclotomic(2040, u).ae(1)
    assert cyclotomic(0,2.5) == 1
    assert cyclotomic(1,2.5) == 2.5-1
    assert cyclotomic(2,2.5) == 2.5+1
    assert cyclotomic(3,2.5) == 2.5**2 + 2.5 + 1
    assert cyclotomic(7,2.5) == 406.234375 

def test_arg_sign():
    assert arg(3) == 0
    assert arg(-3).ae(pi)
    assert arg(j).ae(pi/2)
    assert arg(-j).ae(-pi/2)
    assert arg(0) == 0
    assert isnan(atan2(3,nan))
    assert isnan(atan2(nan,3))
    assert isnan(atan2(0,nan))
    assert isnan(atan2(nan,0))
    assert isnan(atan2(nan,nan))
    assert arg(inf) == 0
    assert arg(-inf).ae(pi)
    assert isnan(arg(nan))
    #assert arg(inf*j).ae(pi/2)
    assert sign(0) == 0
    assert sign(3) == 1
    assert sign(-3) == -1
    assert sign(inf) == 1
    assert sign(-inf) == -1
    assert isnan(sign(nan))
    assert sign(j) == j
    assert sign(-3*j) == -j
    assert sign(1+j).ae((1+j)/sqrt(2)) 

def pairwiseDistances(u, v):
    """
    Pairwise distances between two arrays.

    :param u: first array 
    :type  u: array
    :param v: second array 
    :type  v: array

    :return: array( len(u) x len(v) ) of double
    :rtype: array
    """
    diag1 = N0.diagonal( N0.dot( u, N0.transpose(u) ) )
    diag2 = N0.diagonal( N0.dot( v, N0.transpose(v) ) )
    dist = -N0.dot( v,N0.transpose(u) )\
         -N0.transpose( N0.dot( u, N0.transpose(v) ) )
    dist = N0.transpose( N0.asarray( list(map( lambda column,a:column+a, \
                                        N0.transpose(dist), diag1)) ) )
    return N0.transpose( N0.sqrt( N0.asarray(
        list(map( lambda row,a: row+a, dist, diag2 ) ) ))) 

def cartesianToPolar( xyz ):
    """
    Convert cartesian coordinate array to polar coordinate array: 
    C{ x,y,z -> r, S{theta}, S{phi} }

    :param xyz: array of cartesian coordinates (x, y, z)
    :type  xyz: array

    :return: array of polar coordinates (r, theta, phi)
    :rtype: array
    """
    r = N0.sqrt( N0.sum( xyz**2, 1 ) )
    p = N0.arccos( xyz[:,2] / r )

    ## have to take care of that we end up in the correct quadrant
    t=[]
    for i in range(len(xyz)):
        ## for theta (arctan)
        t += [math.atan2( xyz[i,1], xyz[i,0] )]

    return N0.transpose( N0.concatenate( ([r],[t],[p]) ) ) 

def distanta():
    """Funcția citește conținutul fișierului istoric.tuxy și
    calculează distanța dintre punctul de origine și poziția
    curentă a cursorului.
    """
    cord_x = 0
    cord_y = 0
    fis = open("istoric.tuxy", "r")
    history = fis.read()
    for line in history.splitlines():
        inst = line.split(' ')
        if inst[0] == "STÂNGA":
            cord_x -= int(inst[1])
        if inst[0] == "DREAPTA":
            cord_x += int(inst[1])
        if inst[0] == "SUS":
            cord_y += int(inst[1])
        if inst[0] == "JOS":
            cord_y -= int(inst[1])
    print(cord_x, ' ', cord_y)
    print(cmath.sqrt(cord_x ** 2 + cord_y ** 2)) 

def test_complex_sqrt_accuracy():
    def test_mpc_sqrt(lst):
        for a, b in lst:
            z = mpc(a + j*b)
            assert mpc_ae(sqrt(z*z), z)
            z = mpc(-a + j*b)
            assert mpc_ae(sqrt(z*z), -z)
            z = mpc(a - j*b)
            assert mpc_ae(sqrt(z*z), z)
            z = mpc(-a - j*b)
            assert mpc_ae(sqrt(z*z), -z)
    random.seed(2)
    N = 10
    mp.dps = 30
    dps = mp.dps
    test_mpc_sqrt([(random.uniform(0, 10),random.uniform(0, 10)) for i in range(N)])
    test_mpc_sqrt([(i + 0.1, (i + 0.2)*10**i) for i in range(N)])
    mp.dps = 15 

def gLR(xc, xl):
    if xl == 0:
        return (4*sqrt(xc)*(1 + 9*xc - 9*xc**2 - xc**3 + 6*xc*(1 + xc)*log(xc))).real
    else:
        return (4*sqrt(xc)*(sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        + 10*xc*sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        + xc**2*sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - 5*xl*sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - 5*xc*xl*sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - 2*xl**2*sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - 12*xc*log(2) - 12*xc**2*log(2) + 24*xc*xl*log(2)
        - 12*xl**2*log(2) - 12*(-1 + xc)*xl**2* log(1 - sqrt(xc))
        - 6*xc*log(xc) - 6*xc**2*log(xc) + 12*xc*xl*log(xc) - 3*xl**2*log(xc)
        - 3*xc*xl**2*log(xc) - 6*xl**2*log(sqrt(xc) - 2*xc + xc**1.5)
        + 6*xc*xl**2* log(sqrt(xc) - 2*xc + xc**1.5) - 6*xl**2*log(xl)
        + 6*xc*xl**2*log(xl) + 12*xc*log(1 + xc - xl
        - sqrt(1 + (xc - xl)**2 - 2*(xc + xl)))
        + 12*xc**2*log(1 + xc - xl - sqrt(1 + (xc - xl)**2 - 2*(xc + xl)))
        - 24*xc*xl*log(1 + xc - xl - sqrt(1 + (xc - xl)**2 - 2*(xc + xl)))
        + 6*xl**2*log(1 + xc - xl - sqrt(1 + (xc - xl)**2 - 2*(xc + xl)))
        + 6*xc*xl**2* log(1 + xc - xl - sqrt(1 + (xc - xl)**2 - 2*(xc + xl)))
        + 6*xl**2*log(1 + xc**2 - xl + sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - xc*(2 + xl + sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))))
        - 6*xc*xl**2* log(1 + xc**2 - xl + sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl))
        - xc*(2 + xl + sqrt(xc**2 + (-1 + xl)**2 - 2*xc*(1 + xl)))))).real 

def bezierLength(points, beginT=0, endT=1, samples=1024):
    # https://en.wikipedia.org/wiki/Arc_length#Finding_arc_lengths_by_integrating
    vec = [points[1]-points[0], points[2]-points[1], points[3]-points[2]]
    dot = [vec[0]@vec[0], vec[0]@vec[1], vec[0]@vec[2], vec[1]@vec[1], vec[1]@vec[2], vec[2]@vec[2]]
    factors = [
        dot[0],
        4*(dot[1]-dot[0]),
        6*dot[0]+4*dot[3]+2*dot[2]-12*dot[1],
        12*dot[1]+4*(dot[4]-dot[0]-dot[2])-8*dot[3],
        dot[0]+dot[5]+2*dot[2]+4*(dot[3]-dot[1]-dot[4])
    ]
    # https://en.wikipedia.org/wiki/Trapezoidal_rule
    length = 0
    prev_value = math.sqrt(factors[4]+factors[3]+factors[2]+factors[1]+factors[0])
    for index in range(0, samples+1):
        t = beginT+(endT-beginT)*index/samples
        # value = math.sqrt(factors[4]*(t**4)+factors[3]*(t**3)+factors[2]*(t**2)+factors[1]*t+factors[0])
        value = math.sqrt((((factors[4]*t+factors[3])*t+factors[2])*t+factors[1])*t+factors[0])
        length += (prev_value+value)*0.5
        prev_value = value
    return length*3/samples

# https://en.wikipedia.org/wiki/Root_of_unity
# cubic_roots_of_unity = [cmath.rect(1, i/3*2*math.pi) for i in range(0, 3)] 

def quadratic_roots(a: int, b: int, c: int) -> Tuple[complex, complex]:
    """
    Given the numerical coefficients a, b and c,
    calculates the roots for any quadratic equation of the form ax^2 + bx + c

    >>> quadratic_roots(a=1, b=3, c=-4)
    (1.0, -4.0)
    >>> quadratic_roots(5, 6, 1)
    (-0.2, -1.0)
    >>> quadratic_roots(1, -6, 25)
    ((3+4j), (3-4j))
    """

    if a == 0:
        raise ValueError("Coefficient 'a' must not be zero.")
    delta = b * b - 4 * a * c

    root_1 = (-b + sqrt(delta)) / (2 * a)
    root_2 = (-b - sqrt(delta)) / (2 * a)

    return (
        root_1.real if not root_1.imag else root_1,
        root_2.real if not root_2.imag else root_2,
    ) 

def compare_BR(wc_wilson, l1, l2):
    scale = flavio.config['renormalization scale'][l1+'decays']
    ll = wcxf_sector_names[l1, l2]
    par = flavio.default_parameters.get_central_all()
    wc_obj = flavio.WilsonCoefficients.from_wilson(wc_wilson, par)
    par = flavio.parameters.default_parameters.get_central_all()
    wc = wc_obj.get_wc(ll, scale, par, nf_out=4)
    alpha = flavio.physics.running.running.get_alpha_e(par, scale, nf_out=4)
    e = sqrt(4 * pi * alpha)
    ml = par['m_' + l1]
    # cf. (18) of hep-ph/0404211
    pre = 48 * pi**3 * alpha / par['GF']**2
    DL = 2 / (e * ml) * wc['Cgamma_' + l1 + l2]
    DR = 2 / (e * ml) * wc['Cgamma_' + l2 + l1].conjugate()
    if l1 == 'tau':
        BR_SL = par['BR(tau->{}nunu)'.format(l2)]
    else:
        BR_SL = 1  # BR(mu->enunu) = 1
    return pre * (abs(DL)**2 + abs(DR)**2) * BR_SL 

def add_globals(self):
    "Add some Scheme standard procedures."
    import math, cmath, operator as op
    from functools import reduce
    self.update(vars(math))
    self.update(vars(cmath))
    self.update({
     '+':op.add, '-':op.sub, '*':op.mul, '/':op.itruediv, 'níl':op.not_, 'agus':op.and_,
     '>':op.gt, '<':op.lt, '>=':op.ge, '<=':op.le, '=':op.eq, 'mod':op.mod, 
     'frmh':cmath.sqrt, 'dearbhluach':abs, 'uas':max, 'íos':min,
     'cothrom_le?':op.eq, 'ionann?':op.is_, 'fad':len, 'cons':cons,
     'ceann':lambda x:x[0], 'tóin':lambda x:x[1:], 'iarcheangail':op.add,  
     'liosta':lambda *x:list(x), 'liosta?': lambda x:isa(x,list),
     'folamh?':lambda x: x == [], 'adamh?':lambda x: not((isa(x, list)) or (x == None)),
     'boole?':lambda x: isa(x, bool), 'scag':lambda f, x: list(filter(f, x)),
     'cuir_le':lambda proc,l: proc(*l), 'mapáil':lambda p, x: list(map(p, x)), 
     'lódáil':lambda fn: load(fn), 'léigh':lambda f: f.read(),
     'oscail_comhad_ionchuir':open,'dún_comhad_ionchuir':lambda p: p.file.close(), 
     'oscail_comhad_aschur':lambda f:open(f,'w'), 'dún_comhad_aschur':lambda p: p.close(),
     'dac?':lambda x:x is eof_object, 'luacháil':lambda x: evaluate(x),
     'scríobh':lambda x,port=sys.stdout:port.write(to_string(x) + '\n')})
    return self 

def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 

def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 

def test_sqrt_real_positive(self):
        self.assertAlmostEqual(complex(2, 1),
                               cmath.sqrt(complex(3, 4))) 

def test_sqrt_imaginary_zero(self):
        self.assertAlmostEqual(complex(0.0, 1.73205),
                               cmath.sqrt(complex(-3, 0))) 

def test_sqrt_real_negative(self):
        self.assertAlmostEqual(complex(1, 2),
                               cmath.sqrt(complex(-3, 4))) 

def test_sqrt_imaginary_negative(self):
        self.assertAlmostEqual(complex(1.0, -2.0),
                               cmath.sqrt(complex(-3, -4))) 

def test_sqrt_real_negative(self):
        self.assertAlmostEqual(complex(1, 2),
                               cmath.sqrt(complex(-3, 4))) 

def test_sqrt_real_zero(self):
        self.assertAlmostEqual(complex(1.41421, 1.41421),
                               cmath.sqrt(complex(0, 4))) 

def test_sqrt_real_positive(self):
        self.assertAlmostEqual(complex(2, 1),
                               cmath.sqrt(complex(3, 4))) 

def quadratic(a: float, b: float, c: float) -> Tuple[complex, complex]:
    ''' Compute the roots of the quadratic equation:

            ax^2 + bx + c = 0

        Written in Python as:

            a*x**2 + b*x + c == 0.0

        For example:

            >>> x1, x2 = quadratic(a=8, b=22, c=15)
            >>> x1
            (-1.25+0j)
            >>> x2
            (-1.5+0j)
            >>> 8*x1**2 + 22*x1 + 15
            0j
            >>> 8*x2**2 + 22*x2 + 15
            0j

    '''
    discriminant = cmath.sqrt(b**2.0 - 4.0*a*c)
    x1 = (-b + discriminant) / (2.0 * a)
    x2 = (-b - discriminant) / (2.0 * a)
    return x1, x2 

def bezierRoots(dists, tollerance=0.0001):
    # https://en.wikipedia.org/wiki/Cubic_function
    # y(t) = a*t^3 +b*t^2 +c*t^1 +d*t^0
    a = 3*(dists[1]-dists[2])+dists[3]-dists[0]
    b = 3*(dists[0]-2*dists[1]+dists[2])
    c = 3*(dists[1]-dists[0])
    d = dists[0]
    if abs(a) > tollerance: # Cubic
        E2 = a*c
        E3 = a*a*d
        A = (2*b*b-9*E2)*b+27*E3
        B = b*b-3*E2
        C = ((A+cmath.sqrt(A*A-4*B*B*B))*0.5) ** (1/3)
        roots = []
        for root in cubic_roots_of_unity:
            root *= C
            root = -1/(3*a)*(b+root+B/root)
            if abs(root.imag) < tollerance and root.real > -param_tollerance and root.real < 1.0+param_tollerance:
                roots.append(max(0.0, min(root.real, 1.0)))
        # Remove doubles
        roots.sort()
        for index in range(len(roots)-1, 0, -1):
            if abs(roots[index-1]-roots[index]) < param_tollerance:
                roots.pop(index)
        return roots
    elif abs(b) > tollerance: # Quadratic
        disc = c*c-4*b*d
        if disc < 0:
            return []
        disc = math.sqrt(disc)
        return [(-c-disc)/(2*b), (-c+disc)/(2*b)]
    elif abs(c) > tollerance: # Linear
        root = -d/c
        return [root] if root >= 0.0 and root <= 1.0 else []
    else: # Constant / Parallel
        return [] if abs(d) > tollerance else float('inf') 

def test_sqrt_imaginary_negative(self):
        self.assertAlmostEqual(complex(1.0, -2.0),
                               cmath.sqrt(complex(-3, -4))) 

def test_sqrt_imaginary_zero(self):
        self.assertAlmostEqual(complex(0.0, 1.73205),
                               cmath.sqrt(complex(-3, 0))) 

def solve(self):
    	a = float(self.edit1.text())
    	b = float(self.edit2.text())
    	c = float(self.edit3.text())
        d = (b**2) - (4*a*c)
        root1 = (-b - cmath.sqrt(d)) / (2 * a)
        root2 = (-b + cmath.sqrt(d)) / (2 * a)
        print(root1)
        print(root2)
        root1=str(root1)
        root2=str(root2)
        self.ansx.setText(root1)
        self.ansy.setText(root2)