Complex cube root
Introduction
edit
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
Similarly, are the corresponding components of
and
|
Implementation
edit
Cos φ from cos 3φeditCalculation of cube roots of Wedit
|
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)
editIntroduction
editLet 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
editThis 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
editFour points
edit
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.
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.
|
Point of intersection
edit
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
|
Implementation
editInitialization
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
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
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)
editIntroduction
editLet 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
editsquared:
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:
From From
Simplify Repeat
From |
An Example
editCalculate cube roots of complex number
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: