Complex cube root

Introduction edit

 
Complex number W = complex number w³.
Origin at point  .
  parallel to   axis.
  parallel to   axis.
 
 
 
 
  By cosine triple angle formula:
 
See "3 cube roots of W" in Gallery below.

Let complex numbers   a   b  and   p   q  Let  

When given   aim of this page is to calculate  


In the diagram complex number   where

  •   are the real and imaginary components of   the rectangular components.
  •   are the modulus and phase of   the polar components.

Similarly,   are the corresponding components of  


 

 

 

 


There are 2 significant calculations:

  and

 

Implementation edit

Cos φ from cos 3φ edit

 
Graph of  
Formula and/or graph are used to calculate   if   is given.

The cosine triple angle formula is:   This formula, of form   permits   to be calculated if   is known.

If   is known and the value of   is desired, this identity becomes:     is the solution of this cubic equation.

In fact this equation has three solutions, the other two being  

 


It is sufficient to calculate only   from   accomplished by the following code:

# python code

cosAfrom_cos3Adebug = 0

def cosAfrom_cos3A(cos3A) :
    cos3A = Decimal(str(cos3A))+0
    if 1 >= cos3A >= -1 : pass
    else :
        print ('cosAfrom_cos3A(cos3A) : cos3A not in valid range.')
        return None
    '''
    if cos3A == 0 :
        A = 90 and cos3A = cosA
    if cos3A == 1 :
        A = 0 and cos3A = cosA
    if cos3A == -1 :
        A = 180 and cos3A = cosA
'''
    if cos3A in (0,1,-1) : return cos3A

    # From the cosine triple angle formula:
    a,b,c,d = 4,0,-3,-cos3A

    # prepare for Newton's method.
    if d < 0 : x = Decimal(1)
    else : x = -Decimal(1)
    count = 31; L1 = [x]; almostZero = Decimal('1e-' + str(prec-2))

    # Newton's method:
    while 1 :
        count -= 1
        if count <= 0 :
            print ('cosAfrom_cos3A(cos3A): count expired.')
            break
        y = a*x*x*x + c*x + d
        if cosAfrom_cos3Adebug :
            print ('cosAfrom_cos3A(cos3A) : x,y =',x,y)
        if abs(y) < almostZero : break
        slope = 12*x*x + c
        delta_x = y/slope
        x -= delta_x
        if x in L1[-1:-5:-1] :
            # This value of x has been used previously.
            print ('cosAfrom_cos3A(cos3A): endless loop detected.')
            break
        L1 += [x]

    return x

When   slope   Within area of interest, maximum absolute value of slope   a rather small value for slope.

Consequently, with only 9 passes through loop, Newton's method produces a result accurate to 200 places of decimals .

There are 3 conditions, any 1 of which terminates the loop:

  • abs(y) very close to 0 (normal termination).
  • count expired.
  • endless loop detected with the same value of x repeated.

Calculation of cube roots of W edit

Calculation of 1 root edit

# python code.

def complexCubeRoot (a,b) :
    '''
p,q = complexCubeRoot (a,b) where:
* a,b are rectangular coordnates of complex number a + bi.
* (p + qi)^3 = a + bi.
'''

    a = Decimal(str(a))
    b = Decimal(str(b))

    if a == b == 0 : return (Decimal(0), Decimal(0))
    if a == 0 : return (Decimal(0), -simpleCubeRoot(b))
    if b == 0 : return (simpleCubeRoot(a), Decimal(0))

    Wmod = (a*a + b*b).sqrt()
    wmod = simpleCubeRoot (Wmod)
    cosWφ = a/Wmod
    coswφ = cosAfrom_cos3A(cosWφ)
    p = coswφ * wmod # wreal
    q = (wmod*wmod - p*p).sqrt() # wimag
    # Resolve ambiguity of q, + or -.
    v1 = - 3*p*p*q + q*q*q
    qp = abs(b + v1 )
    qn = abs(b - v1)
    if qp > qn : q *= -1

    return p,q

For function simpleCubeRoot() see Cube_root#Implementation

Calculation of 3 roots edit

See Cube roots of unity.

The cube roots of unity are :  

When   is known, the other 2 cube roots are:

  •  
  •  
# python code
def complexCubeRoots (a,b) :
    '''
r0,r1,r2 = complexCubeRoots (a,b) where:
* a,b are rectangular coordnates of complex number a + bi.
* (p + qi)^3 = a + bi.
* r0 = (p0,q0)
* r1 = (p1,q1)
* r2 = (p2,q2)
'''
    p,q = complexCubeRoot (a,b)
    r3 = Decimal(3).sqrt() # Square root of 3.
    pr3,qr3 = p*r3, q*r3
#  r1 = ((-p-q*r)/2, (p*r - q)/2)
#  r2 = ((-p+q*r)/2, (-p*r - q)/2)
    r0 = (p,q)
    r1 = ((-p-qr3)/2, (pr3 - q)/2)
    r2 = ((-p+qr3)/2, (-pr3 - q)/2)
    return r0,r1,r2

Testing results edit

 

 

# Python code used to test results of code above,

def complexCubeRootsTest (a,b) :
    print ('\n++++++++++++++++++++')
    print ('a,b =', a,b)
    almostZero = Decimal('1e-' + str(prec-5))
    r0,r1,r2 = complexCubeRoots (a,b)
    for root in (r0,r1,r2) :
        p,q = root
        print ('  pq =',(p), (q))
        a_,b_ = (p*p*p - 3*p*q*q), (3*p*p*q - q*q*q)
        if a :
            v1 = abs((a_-a)/a)
            if v1 > almostZero : print ('error *',a_,a,v1)
        else :
            v1 = abs(a_)
            if v1 > almostZero : print ('error !',a_,a,v1)
        if b :
            v1 = abs((b_-b)/b)
            if v1 > almostZero : print ('error &',b_,b,v1)
        else :
            v1 = abs(b_)
            if v1 > almostZero : print ('error %',b_,b,v1)
    return r0,r1,r2
    
import decimal
Decimal = D = decimal.Decimal
prec = decimal.getcontext().prec # precision

cosAfrom_cos3Adebug = 1

for p in range (-10,11,1) :
    for q in range (-10,11,1) :
        a = p*p*p - 3*p*q*q
        b = 3*p*p*q - q*q*q
        complexCubeRootsTest(a,b)

Gallery edit

Method #2 (Graphical) edit

Introduction edit

 
Graphs of 2 curves showing complex cube roots as points of intersection of the 2 curves.

Let complex number  

Then  

Let  

Then:

  and

 


When   is given and   is desired,   may be calculated from the solutions of 2 simultaneous equations:

  and

 


For example, let  


Then equations   and   become (for graphical purposes):

  black curve in diagram, and

  red curve in diagram.


Three points of intersection of red and black curves are:

 

  and

 

interpreted as the three complex cube roots of   namely:

 

  and

 

Proof:

# python code
# Each cube root cubed.
>>> w0 = (-18 + 29j)
>>> w1 = (34.11473670974872 + 1.0884572681198943j)
>>> w2 = (-16.11473670974872 - 30.088457268119896j)
>>> for w in (w0,w1,w2) : w**3
... 
(39582+3799j)
(39582+3799j)
(39582+3799j)
# The moduli of all 3 cube roots.
>>> for w in (w0,w1,w2) : (w.real**2 + w.imag**2) ** 0.5
... 
34.132096331752024
34.132096331752024
34.132096331752024

Preparation edit

This method depends upon selection of most appropriate quadrant.

In the example above, quadrant   is chosen because any non-zero positive value of   intersects red curve and any non-zero negative value of   intersects black curve.

Figures 1-4 below show all possibilities of   and  

Description of method edit

Four points edit

 
Graphs of 2 curves showing 4 points that enclose one of the complex cube roots.
Area enclosed by the 4 points becomes progressively smaller and smaller until the point of intersection is identified.

Assume that   in which case both   are negative and quadrant   is chosen.

In quadrant   any non-zero positive value of x intersects black curve and any non-zero negative value of y intersects red curve.


Choose any convenient negative, non-zero value of  

Let  

Using this value of   calculate coordinates of point   on red curve.

Using   coordinate of point   calculate coordinates of point   on black curve.

Using   coordinate of point   calculate coordinates of point   on red curve.

Using   coordinate of point   calculate coordinates of point   on black curve.


Points   enclose the point of intersection of the 2 curves.

Point of intersection edit

 
Graphs of 2 curves showing complex cube root enclosed by four points  .
Point   intersection of lines   is close to complex cube root, and is starting point for next iteration.

Calculate equations of lines  

Calculate coordinates of point   intersection of lines  


Point   is used as starting point for next iteration.


Area of quadrilateral   becomes smaller and smaller until complex cube root, intersection of red and black curves, is identified.

Implementation edit

Initialization edit

# python code

import decimal
D = decimal.Decimal

precision = decimal.getcontext().prec = 100

useDecimal = 1

Tolerance = 1e-14
if useDecimal :
    Tolerance = D('1e-' + str(decimal.getcontext().prec - 2))


def line_mc (point1,point2) :
    '''
m,c = line_mc (point1,point2)
where y = mx + c.
'''
    x1,y1 = point1
    x2,y2 = point2
    m = (y2-y1) / (x2-x1)
    # y = mx + c
    c = y1 - m*x1
    return m,c


def intersectionOf2Lines (line1, line2) :
    m1,c1 = line1
    m2,c2 = line2
    # y = m1x + c1
    # y = m2x + c2
    # m1x + c1 = m2x + c2
    # m1x - m2x = c2 - c1
    x = (c2-c1)/(m1-m2)
    y = m1*x + c1
    return x,y

def almostEqual (v1,v2,tolerance = Tolerance) :
    '''
status = almostEqual (v1,v2)

For floats, tolerance is 1e-14.

1234567.8901234567 and
1234567.8901234893
are not almostEqual.

1234567.8901234567 and
1234567.8901234593
are almostEqual.
'''
    return abs(v1-v2) < tolerance*abs((v1+v2)/2)

Function two_points() edit

 
Graphs of 2 curves showing 4 points that enclose one of the complex cube roots.
On first invocation function two_points() returns points  
On second invocation function two_points() returns points  
If distance   or distance   is very small, function two_points() returns status True.
def two_points (y,a,b,quadrant) :
    '''
[point1,point2],status = two_points (y,a,b)
'''

    L1 = [] ; yInput = y
    if 0 :
        print ('two_points():')
        s1 = '  a,b,y' ; print (s1, eval(s1))

# a = ppp - pqq - 2pqq = ppp - 3pqq (1)
# b = 2ppq + ppq - qqq = 3ppq - qqq (2)
#
# Using (2)
# 3xxy - yyy = b
#
#     yyy + b
# X = ----------
#        3y
    X = (y**3 + b) / (3*y)

    if isinstance(X,D) : x = X.sqrt()
    else : x = X ** 0.5
    if quadrant in (2,3) : x *= -1
    L1 += [(x,y)]

# Using (1)
# xxx - 3xyy = a
#
#     xxx - a
# Y = -----------
#        3x
#
    Y = (x**3 - a) / (3*x)

    if isinstance(Y,D) : y = Y.sqrt()
    else : y = Y ** 0.5
    if quadrant in (3,4) : y *= -1

    L1 += [(x,y)]

    return L1, almostEqual(y, yInput)

Function pointOfIntersection() edit

 
Graphs of 2 curves showing complex cube root enclosed by four points  .
From points   function pointOfIntersection() calculates coordinates of point  
If distance   is very small, point   is returned as equivalent to intersection of red and black curves.
If distance   is very small, point   is returned as equivalent to intersection of red and black curves.
def pointOfIntersection(y,a,b,quadrant) :
    '''
pointE, status = pointOfIntersection(y,a,b)
y is Y coordinate of point A.
'''
#    print('\npointOfIntersection()')
    t1,status = two_points (y,a,b,quadrant)
    ptA,ptB = t1
    if status :
        # Distance between ptA and ptB is very small.
        # ptA is considered equivalent to the complex cube root.
        return ptA,status

    t2,status = two_points (ptB[1],a,b,quadrant)
    ptC,ptD = t2
    if status :
        # Distance between ptC and ptD is very small.
        # ptC is considered equivalent to the complex cube root.
        return ptC,status

    lineAC = line_mc (ptA,ptC)
    lineBD = line_mc (ptB,ptD)
    pointE = intersectionOf2Lines (lineAC, lineBD)
    return pointE,False

Execution edit

def CheckMake_pq(a,b,x,y) :
    '''
p,q = CheckMake_pq(a,b,x,y)
p = x and q = y.
p,q are checked as valid within tolerance, and are reformatted slightly to improve appearance.
'''
#    print ('\nCheckMake_pq(a,b,x,y)')
#    s1 = '(a,b)' ; print (s1,eval(s1))
#    s1 = '   x'  ; print (s1,eval(s1))
#    s1 = '   y'  ; print (s1,eval(s1))
    P = x**2 ; Q = y**2 ; p=q=-1
    for p in (x,) :
        # a = ppp - 3pqq; a + 3pqq should equal ppp.
        if not almostEqual (P*p, a + 3*p*Q, Tolerance*100) : continue
        for q in  (y,) :
            # b = 3ppq - qqq; b + qqq should equal 3ppq
            if not almostEqual (3*P*q, b + q*Q) : continue
            # Following 2 lines improve appearance of p,q
            if useDecimal :
                # 293.00000000000000000000000000000000000000034 becomes 293
                p,q = [ decimal.Context(prec=precision-3).create_decimal(s).normalize()
                        for s in (p,q) ]
            else :
                # 123.99999999999923 becomes 124.0
                p,q = [ float(decimal.Context(prec=14).create_decimal(s)) for s in (p,q) ]
            return p,q
    # If code gets to here there is internal error.
    s1 = '   p'  ; print (s1,eval(s1))
    s1 = '   q'  ; print (s1,eval(s1))
    s1 = 'P*p, a + 3*p*Q' ; print (s1,eval(s1))
    s1 = '3*P*q, b + q*Q' ; print (s1,eval(s1))
    1/0


def ComplexCubeRoot (a,b, y = 100, count_ = 20) :
    '''
p,q = ComplexCubeRoot (a,b, y, count_)
(p+qi)**3 = (a+bi)
'''
    print ('\nComplexCubeRoot(): a,b,y,count_ =',a,b,y,count_)

    if useDecimal : y,a,b = [ D(str(v)) for v in (y,a,b) ]

    if a == 0 :
        if b == 0 : return 0,0
        if useDecimal : cubeRoot = abs(b) ** (D(1)/3)
        else : cubeRoot = abs(b) ** (1/3)
        if b > 0 : return 0,-cubeRoot
        return 0,cubeRoot
    if b == 0 :
        if useDecimal : cubeRoot = abs(a) ** (D(1)/3)
        else : cubeRoot = abs(a) ** (1/3)
        if a > 0 : return cubeRoot,0
        return -cubeRoot,0

    # Select most appropriate quadrant.
    if a > 0: setx = {2,3}
    else: setx = {1,4}

    if b > 0: sety = {1,2}
    else: sety = {3,4}

    quadrant, = setx & sety
    # Make sign of y appropriate for this quadrant.
    if quadrant in (1,2) : y = abs(y)
    else : y = -abs(y)

    s1 = '    quadrant' ;    print (s1, eval(s1))

    for count in range (0,count_) :
        pointE,status = pointOfIntersection(y,a,b,quadrant)
        s1 = 'count,status' ; print (s1, eval(s1))
        s1 = '    pointE[0]' ;    print (s1, eval(s1))
        s1 = '    pointE[1]' ;    print (s1, eval(s1))
        x,y = pointE
        if status : break

    p,q = CheckMake_pq(a,b,x,y)

    return p,q

An Example edit

p,q = 18,-29
w0 = p + q*1j
W = w0**3
a,b = W.real, W.imag
s1 = '\na,b' ; print (s1, eval(s1))
print ('Calculate one cube root of W =', W)

p,q = ComplexCubeRoot (a, b, -18)
s1 = '\np,q' ; print (s1, eval(s1))

sign = ' + '
if q < 0 : sign = ' - '
print ('w0 = ', str(p) ,sign, str(abs(q)),'i',sep='')
a,b (-39582.0, -3799.0)
Calculate one cube root of W = (-39582-3799j)

ComplexCubeRoot(): a,b,y,count_ = -39582.0 -3799.0 -18 20
    quadrant 4
count,status (0, False)
    pointE[0] 18.29530595866769796981147594954794157453427770441979517949705190002312860122683512802090262517713985
    pointE[1] -29.17605851188829785804056660826030733025475591125914094664311767722817017040959484906193571713185723
count,status (1, False)
    pointE[0] 18.00005338608833140244623119091867294731079673210031643698475740639013316464725974710029414039192178
    pointE[1] -29.00004871113589733281025047965490410760310240487028781197782902118493613318468122514940700337153822
count,status (2, False)
    pointE[0] 18.00000000000411724901243622639913339568271402799883998577879934397861645683113192688877413853607310
    pointE[1] -29.00000000000374389488957142528693701977412078931714190614174861872117809632376941851192785888688849
count,status (3, False)
    pointE[0] 18.00000000000000000000000002432197332441306371168312765664226520285586754485200564671348858119116700
    pointE[1] -29.00000000000000000000000002211642401405406173668734741733917293863102998696594080783163691859719835
count,status (4, False)
    pointE[0] 18.00000000000000000000000000000000000000000000000000000084875301166228708732099292145747224660296019
    pointE[1] -29.00000000000000000000000000000000000000000000000000000077178694502906372553368739653233322951192943
count,status (5, False)
    pointE[0] 18.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
    pointE[1] -29.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
count,status (6, True)
    pointE[0] 18
    pointE[1] -29.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

p,q (Decimal('18'), Decimal('-29'))
w0 = 18 - 29i

Method #3 (Algebraic) edit

Introduction edit

Let complex number  

Then  

Let  

Then:

  and

 


When   is given and   is desired,   may be calculated from the solutions of 2 simultaneous equations   and   where   are known values, and   are desired.

Implementation edit

  squared:  

From  


 

 

 

 

 

 

Let  

 


For   in   substitute  

 


Expand   simplify, gather like terms and result is:

  where:

 

 

 

 

 


Calculate one real root of  

 

From  


Check   against   to resolve ambiguity of sign of  


  is one cube root of  

  may be calculated without ambiguity as follows:


  and

 

From  

From  


Let:

 

 

 

 


Then   become:

 

 

 

 

 

Simplify  

Repeat  

 

 

 

From  

An Example edit

Calculate cube roots of complex number  

 
Graph of   shown as graph of   and showing three values of  .
  axis compressed for clarity.
# python code:

a,b = 39582,3799

s = 64
t = 48*b
u = -(15*b**2 + 27*a**2)
v = b**3

s,t,u,v
(64, 182352, -42518323563, 54828691399)

Calculate roots of cubic function:      

Three roots are:  

# python code:

values_of_Q = Q1,Q2,Q3 = -27239.53953801976, 1.2895380197588122, 24389
q1,q2,q3 = [ abs(Q)**(1/3) for Q in values_of_Q ]
q1 *= -1
values_of_q = q1,q2,q3
s1 = 'values_of_q' ; print(s1, eval(s1))

for m in zip(values_of_q, values_of_Q) :
    q,Q = m
    p = ((b + Q)/(3*q)) ** .5
    sum = p**3 - 3*p*q**2 - a
    if abs(sum) > 1e-10 : p *= -1
    w = p + q*1j
    s1 = 'w, w**3' ; print (s1, eval(s1))
values_of_q (-30.08845726811989, 1.0884572681198943, 29)
w, w**3 ((-16.114736709748723 - 30.08845726811989j  ), (39582+3799j))
w, w**3 (( 34.11473670974874  +  1.0884572681198943j), (39582+3799j))
w, w**3 ((-18.0               + 29.0j               ), (39582+3799j))

Three cube roots of   are: