Objective

Learn about Python's Decimal module.
Learn about Decimal objects.
Review examples of code containing Decimal objects.
Learn how to use Decimal objects to perform calculations on complex numbers.
Become familiar with numbers containing many thousands of decimal digits.
Create a new class CompleX.
Use class CompleX and Python's Decimal module to solve the cubic equation.

Lesson

Calculations using Python's Decimal module can be performed with (almost) any precision selectable by the user. Unfortunately many functions which you might expect to find in the module don't exist, for example, trigonometric functions and functions that manipulate complex numbers. For these functions you have to find them elsewhere or write them for yourself.

Here are some examples that highlight the power of Python's Decimal Module.

Note the following punctuation:

"decimal" is name of module.

"Decimal" is object defined within module "decimal".

Conversion to Decimal object

Several different objects are convertible to Decimal object. Type int converts exactly:

>>> from decimal import *
>>> d1 = Decimal(12345678901123456789011234567890112345678901123456789011234567890112345678901);d1
Decimal('12345678901123456789011234567890112345678901123456789011234567890112345678901')
>>>

Precision is not applied when the Decimal object is initialized or displayed. Default precision of 28 is applied when the Decimal object is used:

>>> d1 + 0; d1 - 0; +d1; -d1
Decimal('1.234567890112345678901123457E+76')
Decimal('1.234567890112345678901123457E+76')
Decimal('1.234567890112345678901123457E+76')
Decimal('-1.234567890112345678901123457E+76')
>>>

Although conversion is not necessary, Decimal object converts exactly.

>>> d2 = Decimal(d1) ; d2
Decimal('123456789011234567890112345678901123456789011234567890112345678901')
>>>

Conversion from float to Decimal is tricky:

>>> f1 = 3.14 ; f1
3.14
>>> d3 = Decimal( f1 ) ; d3 ; +d3
Decimal('3.140000000000000124344978758017532527446746826171875')
Decimal('3.140000000000000124344978758')
>>>

Conversion from float is accurate if correct precision for floats is applied.

>>> getcontext().prec = 14
>>> f1 = 1.13-1.1 ; f1
0.029999999999999805
>>> Decimal(f1); +Decimal(f1); 
Decimal('0.029999999999999804600747665972448885440826416015625')
Decimal('0.030000000000000')
>>>

Also, conversion from float is accurate if float is correctly formatted.

>>> f1 = 1.13-1.1 ; f1
0.029999999999999805
>>> Decimal( '{:1.15f}'.format(f1) )
Decimal('0.030000000000000')
>>>

For simple, accurate conversion to Decimal, convert float to str first:

>>> f1 = 3.14 ; str(f1) ; Decimal(str(f1))
'3.14'
Decimal('3.14')
>>> 
>>> Decimal( '{}'.format(f1) )
Decimal('3.14')
>>>

Conversion from complex to Decimal:

>>> Decimal(3+4j)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: conversion from complex to Decimal is not supported
>>>

Conversion from str representing number to Decimal is accurate.

>>> Decimal('1234.5678901e3')
Decimal('1234567.8901')
>>>

eval() is more forgiving than Decimal().

>>> s1
'  -  3456e-3  '
>>> Decimal(s1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
decimal.InvalidOperation: [<class 'decimal.ConversionSyntax'>]
>>> Decimal( str(eval(s1)) )
Decimal('-3.456')
>>>

Conversion from DecimalTuple is accurate:

>>> dt1 = DecimalTuple(sign=1, digits=(0,0,3, 4, 5, 6,7,8,0,0), exponent=-4) ; dt1 ; Decimal(dt1)
DecimalTuple(sign=1, digits=(0, 0, 3, 4, 5, 6, 7, 8, 0, 0), exponent=-4)
Decimal('-3456.7800') # Leading zeroes are dropped. Trailing zeroes are kept.
>>>

DecimalTuple may be created with bytes for digits:

>>> dt1 = DecimalTuple(sign=1, digits=bytes((0,0,3, 4, 5, 6,7,8,0,0)), exponent=-4) ; dt1
DecimalTuple(sign=1, digits=b'\x00\x00\x03\x04\x05\x06\x07\x08\x00\x00', exponent=-4)

DecimalTuple may not be used with bytes for digits:

>>> Decimal(dt1)
Traceback (most recent call last):
  File "<pyshell#43>", line 1, in <module>
    Decimal(dt1)
ValueError: coefficient must be a tuple of digits

Conversion from well formed list or tuple is accurate:

>>> t1=(1,(2,3,4,5,6),-4);
>>> t2=(1,list(t1[1]),-4);
>>> L1=list(t1);
>>> L2=list(t2);
>>> t1;t2;L1;L2;
(1, (2, 3, 4, 5, 6), -4)
(1, [2, 3, 4, 5, 6], -4)
[1, (2, 3, 4, 5, 6), -4]
[1, [2, 3, 4, 5, 6], -4]
>>> [ Decimal(v) for v in (t1,t2,L1,L2) ]
[Decimal('-2.3456'), Decimal('-2.3456'), Decimal('-2.3456'), Decimal('-2.3456')]
>>> 
>>> { Decimal(v) for v in (t1,t2,L1,L2) }
{Decimal('-2.3456')}
>>>

Using Decimal objects

Decimal objects work well with the usual arithmetic operators:

>>> from decimal import *
>>> D = Decimal
>>> D('234.567') + D('.000000000000000456')
Decimal('234.567000000000000456')
>>> D('234.567') - D('.000000000000000456')
Decimal('234.566999999999999544')
>>> D('234.567') * D('.000000000000000456')
Decimal('1.06962552E-13')
>>> D('234.567') / D('.000000000000000456')
Decimal('514401315789473684.2105263158')
>>> D('234.567') ** ( D('1')/D(3)  )
Decimal('6.167213327076116022863000610') # Cube root.
>>> D('234.567') % 1
Decimal('0.567') # the fractional part.
>>> D('-234.567') % 1
Decimal('-0.567')
>>> D('-45234.567') % 360
Decimal('-234.567')
>>>

If you are doing much heavy math containing cube roots, it might be advantageous for you to write your own cube root function using Newton's method, for example. Newton's method is much faster than raising a number to the power (1/3).

Decimal objects work well with many of the familiar built-in functions:

>>> d1 = D(-5);d1
Decimal('-5')
>>> abs(d1)
Decimal('5')
>>> 
>>> Decimal( ascii(123.456) )
Decimal('123.456')
>>> 
>>> bool(D(4))
True
>>> bool(D(0))
False
>>> 
>>> complex( D('34.56') ) ;  complex(4, D(3))
(34.56+0j)
(4+3j)
>>> 
>>> divmod ( -D('2345678.0987654321'), 360 )
(Decimal('-6515'), Decimal('-278.0987654321'))
>>> divmod ( -D('2345678.0987654321'), 1 )
(Decimal('-2345678'), Decimal('-0.0987654321'))
>>> 
>>> float ( -D('2345678.0987654321'))
-2345678.098765432
>>> 
>>> int ( -D('2345678.987654321'))
-2345678
>>>
>>> isinstance( D(10), Decimal )
True
>>> type( D(10) )
<class 'decimal.Decimal'>
>>> 
>>> max(100, -23, D(44))
100
>>> min(100, -23, D(44))
-23
>>> 
>>> pow(3,D(2))
Decimal('9')
>>> 
>>> sorted(( 3,45,-100, D('234.56') ))
[-100, 3, 45, Decimal('234.56')]
>>> 
>>> str(D('456.78'))
'456.78'
>>> 
>>> sum(( 3,45,-100, D('234.56') ))
Decimal('182.56')
>>>

Decimal objects and attributes

>>> D('3.14159').as_integer_ratio()
(314159, 100000)
>>> D(3.14159).as_integer_ratio()
(3537115888337719, 1125899906842624)
>>>
>>> D('3.14159').as_tuple()
DecimalTuple(sign=0, digits=(3, 1, 4, 1, 5, 9), exponent=-5)
>>> D(3.14159).as_tuple()
DecimalTuple(sign=0, digits=(3, 1, 4, 1, 5, 8, 9, 9, 9, 9, 9, 9, 9, 9, ............... , 7, 0, 9, 9, 6, 0, 9, 3, 7, 5), exponent=-50)
>>> 
>>> D(3.14159).compare(D('3.14159'))
Decimal('-1') # D(3.14159) < D('3.14159')
>>> 
>>> D('-3.14159').copy_abs()
Decimal('3.14159')
>>> abs(D('-3.14159'))
Decimal('3.14159')
>>> 
>>> D('-3.14159').is_normal()
True
>>> D('-3.14159').is_zero()
False
>>> D('-3.14159').is_infinite()
False
>>> 
>>> D(7).max(D(-9))
Decimal('7')
>>> D(7).max_mag(D(-9))
Decimal('-9')
>>> D(7).min(D(-9))
Decimal('-9')
>>> D(7).min_mag(D(-9))
Decimal('7')
>>> 
>>> Decimal('1.41421356').quantize(Decimal('1.000'))
Decimal('1.414')
>>> Decimal('1.41451356').quantize(Decimal('1.000'))
Decimal('1.415')
>>> Decimal('1.41421356').quantize(Decimal('.001'))
Decimal('1.414')
>>> Decimal('1.41451356').quantize(Decimal('.001'))
Decimal('1.415')
>>> 
>>> Decimal('0.321000e+2').normalize()
Decimal('32.1')
>>> Decimal('3.2100e1').normalize()
Decimal('32.1')
>>> Decimal('32100.00000e-3').normalize()
Decimal('32.1')
>>>

Exponential operations

>>> D('1').exp()
Decimal('2.718281828459045235360287471') # Value of e, base of natural logarithms.
>>> D('2').exp()
Decimal('7.389056098930650227230427461') # e ** 2
>>> D('3.14159').ln() # Natural logarithm (base e).
Decimal('1.144729041185178381216412580') # e ** 1.144729... = 3.14159
>>> (D('1234.5678').ln()).exp()
Decimal('1234.567800000000000000000000') # This simple test gives an impressive result.
>>> (D('1234.5678').exp()).ln()
Decimal('1234.567800000000000000000000') # This also.
>>> 
>>> # Raising number to power:
>>> D('1.6') ** D('2.3')
Decimal('2.947650308163391181711649979')
>>> D('1.6').ln()*D('2.3').exp()
Decimal('4.687901952522058518151002058')
>>> (D('1.6').ln()*D('2.3')).exp() # Parentheses are important.
Decimal('2.947650308163391181711649980')
>>> 
>>> D('1.6').sqrt() # Method .sqrt() is very fast.
Decimal('1.264911064067351732799557418')
>>>

Logical operations

The Decimal object D('1110010100011') has the appearance of a binary number and logical operations act on Decimal objects whose only digits are 0,1.

>>> D(1100020011).logical_and(D(1111))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
decimal.InvalidOperation: [<class 'decimal.InvalidOperation'>]
>>>

When using python's decimal module for logical operations on integers, convert int to appropriate Decimal object first.

>>> int1 = 23 ; bin(int1)
'0b10111'
>>> D( bin(int1)[2:] )
Decimal('10111')
>>>

After the logical operation convert Decimal object with appearance of binary number to int.

>>> int(str(D('10111')),2)
23
>>>

>>> D(110001111).logical_invert()
Decimal('1111_1111_1111_1111_1110_0111_0000') # Default precision of 28.
>>> # Same as:
>>> D(110001111).logical_xor(D('1'*getcontext().prec))
Decimal('1111_1111_1111_1111_1110_0111_0000')
>>>

>>> D(110001111).logical_and(D(1111))
Decimal('1111')
>>> # Equivalent to:
>>> bin(int('18F',16) & int('1111',2))
'0b1111'
>>>

>>> D(1_1000_1111).logical_xor(D('1100_1100'))
Decimal('1_0100_0011')
>>>

>>> D(110001111).shift(3) # +ve means shift left.
Decimal('110001111000')
>>> D(110001111).shift(-3) # -ve means shift right.
Decimal('110001') # Bits shifted out to the right are lost.
>>> getcontext().prec=10
>>> D(110001111).shift(3)
Decimal('1111000') # Bits shifted out to the left are lost.
>>>

>>> getcontext().prec=10
>>> D(110001111).rotate(3) # To left.
Decimal('1111011')
>>> D(110001111).rotate(-3) # To right.
Decimal('1110110001')
>>> # Same as:
>>> D('0110001111').logical_and(D(111)).shift(7) + D('0110001111').shift(-3)
Decimal('1110110001')
>>>

Context objects and attributes

Contexts are environments for arithmetic operations. They govern precision, set rules for rounding, determine which signals are treated as exceptions, and limit the range for exponents.

After the decimal module has been imported, three supplied contexts are available:

>>> from decimal import *
>>>
>>> DefaultContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> BasicContext
Context(prec=9, rounding=ROUND_HALF_UP, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[Clamped, InvalidOperation, DivisionByZero, Overflow, Underflow])
>>>
>>> ExtendedContext
Context(prec=9, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[])
>>>

After the decimal module has been imported, the current context is the same as DefaultContext.

>>> from decimal import *
>>>
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> str(getcontext())
'Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])'
>>>
>>> str(getcontext()) == str(DefaultContext)
True
>>>

After importing the decimal module, set the current context (if necessary) as appropriate for your planned use of the decimal module.

>>> from decimal import *
>>> getcontext().prec = 20
>>> getcontext().clear_flags()
>>> getcontext()
Context(prec=20, rounding=ROUND_HALF_UP, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[Clamped, InvalidOperation, Overflow, Underflow])
>>>

For a list of valid signals:

>>> getcontext().flags['']=False
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
KeyError: 'valid values for signals are:\n  [InvalidOperation, FloatOperation, DivisionByZero,\n   Overflow, Underflow, Subnormal, Inexact, Rounded,\n   Clamped]'
>>>

To create a new context copy an existing context:

>>> myContext = BasicContext # This creates a shallow copy.
>>>
>>> myContext = BasicContext.copy()
>>> myContext.prec = 88
>>> myContext ; BasicContext
Context(prec=88, rounding=ROUND_HALF_UP, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[Clamped, InvalidOperation, DivisionByZero, Overflow, Underflow]) # Deep copy.
Context(prec=9, rounding=ROUND_HALF_UP, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[Clamped, InvalidOperation, DivisionByZero, Overflow, Underflow])
>>>

or use the constructor Context():

>>> from decimal import *
>>>
>>> myContext = Context() ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> str(myContext) == str(DefaultContext)
True
>>> myContext = Context(rounding=ROUND_HALF_UP,flags=[]) ; myContext
Context(prec=28, rounding=ROUND_HALF_UP, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> myContext = Context(Emax=9999, flags=[], traps=[InvalidOperation, DivisionByZero]) ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=9999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero])
>>>

To modify a context:

>>> myContext.Emin = -9999 ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-9999, Emax=9999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero])
>>>
>>> myContext.flags[Inexact] = True ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-9999, Emax=9999, capitals=1, clamp=0, flags=[Inexact], traps=[InvalidOperation, DivisionByZero])
>>>
>>> myContext.traps = DefaultContext.traps ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-9999, Emax=9999, capitals=1, clamp=0, flags=[Inexact], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> myContext.clear_flags() ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-9999, Emax=9999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>>
>>> myContext.clear_traps() ; myContext
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-9999, Emax=9999, capitals=1, clamp=0, flags=[], traps=[])
>>>

To set current context:

>>> from decimal import *
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])
>>> setcontext(Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[], traps=[]))
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[], traps=[])
>>>

The rules are:

A bad operation sets a flag.
A good operation does not clear a flag.
Before a decimal operation set traps and clear flags.
After a decimal operation process exceptions and check flags.

>>> getcontext().clear_flags() ; getcontext().clear_traps() ; getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[])
>>>
>>> 1/0
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: division by zero
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[])
>>>
>>> Decimal(1)/0
Decimal('Infinity')
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[DivisionByZero], traps=[])
>>> 
>>> getcontext().clear_flags() ; getcontext().clear_traps()
>>>
>>> 1234567890123456789012345678901234567890+2
1234567890123456789012345678901234567892
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[])
>>>
>>> 1234567890123456789012345678901234567890+Decimal(2)
Decimal('1.234567890123456789012345679E+39')
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[])
>>>

Methods

Many of the methods available for Context objects have the same names as corresponding methods available for Decimal objects, for example : max, quantize, shift and sqrt.

However they usually take an extra argument so that they make sense when attached to Context object;

>>> Decimal(3).sqrt()
Decimal('1.7320508075688772935274463415058723669428053')
>>> BasicContext.sqrt(Decimal(3))
Decimal('1.73205081')
>>>

Others such as clear_traps(), clear_flags() make sense only when attached to Context object.

Context objects are useful if you want to perform an arithmetic operation in a temporary environment without changing current environment.

>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[])
>>> Context(prec=99).ln(Decimal(5))
Decimal('1.60943791243410037460075933322618763952560135426851772191264789147417898770765776463013387809317961')
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[])
>>> 
.......
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[], traps=[])
>>> Context(prec=14).create_decimal_from_float(1.13 - 1.1)
Decimal('0.030000000000000')
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[], traps=[])
>>>

Also if you want to apply a trap to a conversion without affecting current environment:

>>> getcontext().clear_flags()
>>> Decimal('123.456').quantize(Decimal('.01'))
Decimal('123.46')
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[])
>>> getcontext().clear_flags()
>>> Context(traps=[Inexact]).quantize( Decimal('123.456'), Decimal('.01') )
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
decimal.Inexact: [<class 'decimal.Inexact'>]
>>> getcontext()
Context(prec=44, rounding=ROUND_HALF_UP, Emin=-999, Emax=999, capitals=1, clamp=0, flags=[], traps=[])
>>>

To check the result of arithmetic operation in a temporary environment:

>>> myContext=Context(prec=14, flags=[], traps=[], rounding=ROUND_HALF_UP, Emax=99, Emin=-99)
>>> myContext
Context(prec=14, rounding=ROUND_HALF_UP, Emin=-99, Emax=99, capitals=1, clamp=0, flags=[], traps=[])
>>> myContext.quantize( Decimal('123.456'), Decimal('.01') )
Decimal('123.46')
>>> myContext
Context(prec=14, rounding=ROUND_HALF_UP, Emin=-99, Emax=99, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[])
>>> 
>>> myContext.flags[Inexact]
True
>>> myContext.flags[Overflow]
False
>>>

Is Decimal object int?

>>> d1 = Decimal('123.000');d1
Decimal('123.000')
>>> myContext.clear_flags()
>>> myContext.quantize( d1, Decimal(1) )
Decimal('123')
>>> myContext.flags[Inexact]
False # Conversion was exact. d1 is equivalent to int.
>>>
>>> d1 = Decimal('123.007');d1
Decimal('123.007')
>>> myContext.clear_flags()
>>> myContext.quantize( d1, Decimal(1) )
Decimal('123')
>>> myContext.flags[Inexact]
True # Conversion was not exact. d1 is not equivalent to int.
>>>

Making the Decimal object

def print_error (error, x=None, thisName=None) :
    """
Prints error derived from sys.exc_info().
    """
    list1 = error.split(',')
    v1 = ','.join(list1[1:-1])
    if thisName : print (thisName)
    if x != None : print ('  Input =', (str(x))[:60])
    print (' ', list1[0])
    print (' ', v1)
    print (' ', list1[-1])

The following function verifies that we are working with Decimal objects. It provides for creation of Decimal object from DecimalTuple.

import decimal

getcontext = decimal.getcontext

getcontext().prec += 5  # Extra digits for intermediate steps.

dD = D = Decimal = decimal.Decimal
DT = decimal.DecimalTuple

def makeDecimal (x, flag=0) :
    '''
output = makeDecimal (x [, flag])
x is a single object convertible to Decimal object.
returns Decimal object or
        None on error
'''

    thisName = 'makeDecimal (x) :'
    if flag :
        print ()
        print (thisName)
        print ('  Input x =', type(x), (str(x))[:60])
        
    if isinstance(x, dD) : return x
    if isinstance(x, int) : return dD(x)
    if isinstance(x, float) : return dD(str(x))
    if isinstance(x, str) :
        try :
            error = ''
            output = dD(x)
        except : error = str(sys.exc_info())
        if error :
            if flag : print_error (error) 
            return None
        return output
        
    if isinstance(x, (list,tuple,DT)) :
        try :
            error = ''
            v1,v2,v3 = x
            if isinstance(v2, (tuple,list)) : output = dD(x)
            else : output = dD( ( v1, list(v2), v3 ) )
        except : error = str(sys.exc_info())
        if error :
            if flag : print_error (error) 
            return None
        return output

    if flag :
        print ('  Input not recognized.')
    return None

Trigonometric Functions

The following trigonometric functions are sufficient to convert from complex type polar to complex type rectangular and vice-versa.

For the functions $\cos x$ and $\sin x$ see Recipes.

arctan x

We'll try $\arctan x$ first because it's easily understood and implemented, and it can be used to calculate $\pi$ because the result of this expression (as in all theoretical math) is in radians.

$\arctan x=\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{2n+1}}x^{2n+1}=x-{\frac {x^{3}}{3}}+{\frac {x^{5}}{5}}-\cdots {\text{for }}|x|\leq 1,\ x\neq \pm i$

def arctan (tanθ) : # The Python interpreter recognizes international text. Handy.
  # value returned is in radians.
  # Check input:
  x = makeDecimal (tanθ)
  if not isinstance(x, Decimal) :
      print ('arctan: type of input should be Decimal.')
      exit(79)
  
  if not x : return 0 # tan(0) = 0

  if abs(x) > 1 :
      # abs() function is valid for Decimal objects.
      print ('arctan: abs(x) should be <= 1.')
      exit(78)
  if x < 0 :
      print ('arctan: input must be >= 0.')
      exit(77) # Only non-negative values in this implementation.

  sum = x
  x_sq = x*x

  divisor = 3

  last_dividend = x

  multiplier = -1

  getcontext().prec += 2  # extra digits for intermediate steps

  almostZero = Decimal('1e-' + str(getcontext().prec))

  while 1 :
    this_dividend = last_dividend * x_sq
    if this_dividend < almostZero : break
    sum += multiplier * this_dividend / divisor

    last_dividend = this_dividend
    multiplier *= -1
    divisor += 2

  getcontext().prec -= 2  # Restore original precision.
  return sum+0 # Apply original precision to sum.

Calculating π

There are many ways to calculate $\pi$ . For example:

$\pi =6*\arctan(\tan 30)=6*\arctan({\frac {\sqrt {3}}{3}}).$

Based on the following, we will perform six calculations of $\pi$ and compare the results.

$\pi =30*\arctan(\tan 6)$

$\pi =24*\arctan(\tan 7.5)$

$\pi =12*\arctan(\tan 15)$

$\pi =8*\arctan(\tan 22.5)$

$\pi =6*\arctan(\tan 30)$

$\pi =5*\arctan(\tan 36)$

precision = getcontext().prec
getcontext().prec = 502

tan6 = (
    ( (D('10') - D('20').sqrt()).sqrt() + D('3').sqrt() - D('15').sqrt() )
                                       /
                                       2
)
tan7_5 = D(6).sqrt() - D(3).sqrt() + D(2).sqrt() - 2
tan15 = D(2) - D(3).sqrt()
tan22_5 = D(2).sqrt() - 1
tan30 = D('3').sqrt()/3
tan36 = ( D(5) - 2*(D(5).sqrt()) ).sqrt()

values = [
    (tan6,   30),
    (tan7_5, 24),
    (tan15,  12),
    (tan22_5, 8),
    (tan30,   6),
    (tan36,   5)
]

values_of_π = [ π
    for v1 in [ D('1e-' + str(getcontext().prec-3)) ]
    for (tanθ, multiplier) in values
        for π in [
            (multiplier*arctan( tanθ )).quantize( v1, rounding=decimal.ROUND_HALF_UP )
        ]
]

print ('number of calculations =', len(values_of_π))

set2 = set(values_of_π)
len(set2) == 1 or exit(93) # All values in values_of_π must be equal.
π, = list(set2)
str_π = str(π)

L3 = [ str1
    for desired_length in [ 70 ]
    for start in range(0, len(str_π), desired_length)
        for str1 in [ str_π[start:start+desired_length] ]
]

str1 = '"{}"'.format( '" + \n"'.join(L3) )

print (
'''
π = (
{} )
len(str(π)) = {}
isinstance(π, Decimal): {}
'''.format ( str1,
             len(str_π),
             isinstance(π, Decimal)
           )
)

getcontext().prec = precision

number of calculations = 6

π = (
"3.14159265358979323846264338327950288419716939937510582097494459230781" +
"6406286208998628034825342117067982148086513282306647093844609550582231" +
"7253594081284811174502841027019385211055596446229489549303819644288109" +
"7566593344612847564823378678316527120190914564856692346034861045432664" +
"8213393607260249141273724587006606315588174881520920962829254091715364" +
"3678925903600113305305488204665213841469519415116094330572703657595919" +
"5309218611738193261179310511854807446237996274956735188575272489122793" +
"81830119491" )
len(str(π)) = 501
isinstance(π, Decimal): True

All six calculations agree to 500 digits of precision and π is available globally as a Decimal object.

tan (Θ / 2)

$\tan {\frac {\theta }{2}}$ $=\csc \theta -\cot \theta$ .

def tanθ_2 (tanθ):
    '''
tan(θ/2) = csc θ - cot θ
x = tanθ
'''
    x = makeDecimal(tanθ)
    
    isinstance(x, Decimal) or ({}['tanθ_2: x should be Decimal.'])

    if x == 0 : return 0
    
    (x < 0) and ({}['tanθ_2: input < 0.'])

    getcontext().prec += 2  # extra digits for intermediate steps
    
    # cscθ = ((Decimal(1)+x*x).sqrt()) / x
    # cotθ = Decimal(1)/x
    # tan_θ2 = cscθ - cotθ
    tan_θ2 = ((1+x*x).sqrt() - 1) / x
    getcontext().prec -=2

    return (tan_θ2 + 0)

decATAN2 (y, x)

For the corresponding function using floating point arithmetic see math.atan2(y, x)

This function invokes $arctan(x)\ 0\leq x\leq 1$ .

If $\tan \theta >1,\ x={\frac {1}{\tan \theta }}$ and the result is ${\frac {\pi }{2}}-\arctan(x).$

However this function invokes $\arctan(x)$ with a value less than $0.325$ so that the expression $x^{2n+1}$ in the expansion of $\arctan x$ vanishes as rapidly as possible. Therefore:

$75=90-\arctan({\frac {1}{\tan 75}})$

$68=4\arctan(\tan {\frac {68}{4}})$

$32=2\arctan(\tan {\frac {32}{2}})$

$12=\arctan(\tan 12)$

def decATAN2 (y,x) :
  '''
               y
Input is value -.
               x
Both x,y must be convertible to Decimal object.
Only 1 of x or y may be 0.
Returns value of angle in degrees.
Value of π must be available globally.
'''
  x = makeDecimal(x)
  if not isinstance(x, Decimal) :
      print ('decATAN2: type of x',type(x),'should be Decimal.')
      exit(70)

  y = makeDecimal(y)
  if not isinstance(y, Decimal) :
      print ('decATAN2: type of y',type(y),'should be Decimal.')
      exit(69)

  if x == y == 0 :
      print ('decATAN2: both x and y cannot be 0.')
      exit(68)

  if y == 0 :
    if x > 0 : return 0
    return 180

  if x == 0 :
    if y > 0 : return 90
    return 270

  if abs(x) == abs(y) :
      if x > 0 :
          if y > 0 : return 45
          return 360-45
      if y > 0 : return 180-45
      return 180+45

  getcontext().prec += 2  # Extra digits for intermediate steps.
  tanθ = abs(y)/abs(x)

  flip = 0
  if tanθ > Decimal('3.078') : # > 72 degrees
    tanθ = 1/tanθ
    flip += 1

  reductionCount = 0
  while tanθ > Decimal('0.325') : # > 18 degrees
    tanθ = tanθ_2 (tanθ)
    reductionCount += 1

  θ = arctan(tanθ)
  if flip :
    θ = π/2 - θ
  else :
    while reductionCount :
      θ *= 2
      reductionCount -= 1

  θinDegrees = θ*180/π
  if x > 0 :
    if y < 0 :
      θinDegrees = 360-θinDegrees
  else :
    if y > 0 :
      θinDegrees = 180-θinDegrees
    else :
      θinDegrees = 180+θinDegrees

  getcontext().prec -= 2
  return θinDegrees+0

degreesToRadians (θinDegrees)

def degreesToRadians (θinDegrees) :
  '''
Value of π must be available globally.
Value returned: -π < radians <= π
-180 degrees is returned as π.
270 degrees is returned as -π/2
'''
  x = makeDecimal(θinDegrees)

  if not isinstance(x, Decimal) :
      print ('degreesToRadians: type of x should be Decimal.')
      exit(54)

  x = x % 360
  if x < 0 : x += 360
  if x > 180 : x -= 360
  return x * π/180

Complex Functions

Within the context of this page a complex number is contained in a tuple with three members, thus:

(modulus, phaseInDegrees, 'polar')

The above tuple represents complex number $modulus\ *\ (\ \cos(phaseInDegrees)\ +\ 1j*\sin(phaseInDegrees)\ )$ .

Or:

(realPart, imaginaryPart, 'rect')

where $realPart,imaginaryPart$ are the rectangular coordinates of complex number $realPart\ +\ 1j*imaginaryPart$ .

The four values $modulus,phaseInDegrees,realPart,imaginaryPart$ are all Decimal objects.

The rectangular format is useful for addition and subtraction of complex numbers.

The polar format is useful for raising a complex number to a power including a power less than unity.

Both formats are useful for multiplication, division and square root.

When working with polar format it's generally more advantageous to work with a positive modulus. Therefore:

$-n*(\cos \phi +1j*\sin \phi )=n*(\cos(\phi \pm 180)+1j*\sin(\phi \pm 180))$

and, for example:

${\sqrt {-4}}={\sqrt {-4(\cos 0+1j*\sin 0)}}={\sqrt {4(\cos 180+1j*\sin 180)}}=2(\cos 90+1j*\sin 90)=2*1j*1=2j$

and the other sqrt of $-4$ is:

${\sqrt {4(\cos(180-360)+1j*\sin(180-360))}}={\sqrt {4(\cos(-180)+1j*\sin(-180))}}=2(\cos(-90)+1j*\sin(-90))=2*1j*(-1)=-2j$ .

>>> (2j)**2 ; (-2j)**2
(-4+0j)
(-4+0j)
>>>

Both of the following complex tuples equate to 0:

>>> (0,0,'rect')
>>> (0,anyNumber,'polar')

The following functions will enable us to do some serious math with complex numbers, such as solving the cubic equation with three real roots.

Administrative functions

check_complex(x)

This function verifies that the object is a valid complex tuple.

def check_complex(x, flag = 0) :
    """
status = check_complex(x [, flag]) 
x must be :
    (v1,v2,'rectangular') or
    (v1,v2,'polar')
    v1,v2 must be type Decimal.
    """
    thisName = 'check_complex (x) :'

    if flag :
        print ()
        print (thisName)
        print ('  Input x =', type(x), (str(x))[:60])

    try:
        error = ''
        v1,v2,str1 = x
        isinstance (v1,dD) or ({}['v1 must be type Decimal.'])
        isinstance (v2,dD) or ({}['v2 must be type Decimal.'])
        (isinstance (str1,str) and str1) or ({}['str1 must be valid string.'])
        str1 = str1.lower()
        if str1 == 'rectangular'[:len(str1)] : pass
        elif str1 == 'polar'[:len(str1)] : pass
        else : ({}['str1 must be "rectangular" or "polar".'])
    except : error = str(sys.exc_info())
    if error :
        if flag : print_error (error)
        return False
    return True

str_to_complex(x)

This function accepts a string like $(2345678901234567.8903456789-9876543219654321098.765432765j)$ and preserves precision of both real and imaginary parts.

import re

# The re pattern for numbers to be found in string of type complex.
# In string like '  1234567890987654323456E-6 + 987654321234567890987e+2J   '
# extracts '1234567890987654323456E-6' and '+987654321234567890987e+2J'.
# eval allows '  12_3456_7890_9876_5432_3456E-6 + 987654321234567890987e+2J   '.
# Conversion to Decimal allows: '__987_65432123____4567890_987___'.
# eval does not allow ' 0123 ' but allows '  2e04 '.
# Conversion to Decimal allows: ' 0123  '
digits  = '[_0123456789]'                        # '_' is included.
integer = '[_0123456789]{1,}'                    # Matches 1 23 345 4567 56_789
exponent = '[Ee][\+\-]{{0,1}}{}'.format(integer) # Matches e8 e+8 e-8 or E8 E+8 E-8
Exponent = '({}){{0,1}}'.format(exponent)        # exponent 0 or 1 times.
float1 = '({})\.{{0,1}}'.format(integer)         # Matches 123  2345.
float2 = '{}{{0,}}\.{}'.format(digits, integer)  # Matches   123.345  .345
float_ = '(({})|({}))'.format(float2,float1)     # float2 before float1.
rvalue = '{}({})'.format(float_,Exponent)        # real value
rvalue_signed = '[\+\-]{{0,1}}{}'.format(rvalue)
ivalue_signed = '{}[Jj]'.format(rvalue_signed)

def str_to_complex (input_string) :
    """
a,b = str_to_complex (input_string)
input_string could be '   1234567890987654323456E-6 + 987654321234567890987e+2J    '
input_string could be ' ' '  ( 12345678909876543234561234567890987.654323456E-6 +
                             987654321234567890987987654321.23764567890987e+2J  )  ' ' '
or ' ' '
    - (
 (
    12345678_90543123456789054_312345678905431234567890543e-014 +
    543123456789054_3123456789054312345_67890543654327453E-1_3J
) )

    ' ' '
or ' ' '
    - (
 (
    12345678_90543123456789054_312345678905431234567890543e-014 + # Note the '_'
    543123456789054_3123456789054312345_67890543654327453E-1_3J
) )

    ' ' '
or '4', '-4', '(4)', '(-4)', '-(4)', '-(-4)'
or '( +3j)', '+(-3J)', '-(3j)', '-(-3J)'
or '( 4+3j)', '+(-3 + 4J)','-( 4+3j)', '-(-3 + 4J)',
This function retains precision of a,b
    """
    
    try :
        status = 0
        isinstance(input_string,str) or ({}['input_string not type str.'])
        str1 = input_string.strip()
        cx1 = eval(str1)
        isinstance(cx1,(int,float,complex)) or ({}['cx1 not desired type.'])
    except : status = 1
    if status : return

    # This code removes white lines and comments, if any, at end of each line.
    new_line = '''
'''[-1:]
    lines = [ line for Line in str1.split(new_line) for line in [ Line.rstrip() ] if line ]
    lines = [ v for line in lines
                for parts in [ line.split('#') ]
                for v in [ parts[0]] ]
    str1 = ''.join(lines)
    str1 = ''.join(str1.split())

    if isinstance (cx1, (int,float)) :
        resultr = re.search(rvalue_signed, str1)
        resultr or ({}['rvalue not found.'])
        dD1 = dD(resultr[0]) ; v1 = eval(str(dD1))
        if v1 == cx1 : return dD1,dD(0)
        if v1 == -cx1 : return -dD1,dD(0)
        ({}['dD1 not recognized.'])

    # cx1 must be complex. It must contain imaginary value.
    resulti = re.search(ivalue_signed, str1)
    resulti or ({}['ivalue not found.'])
    str2j = resulti[0] ; str2 = str2j[:-1]
    dD2 = dD(str2)
    str1 = ''.join(str1.split(str2j))

    # cx1 may contain real value.
    resultr = re.search(rvalue_signed, str1)
    if resultr : dD1 = dD(resultr[0])
    else : dD1 = dD(0)

    cx2 = complex(dD1,dD2)
    if cx2.real == cx1.real : pass
    elif cx2.real == -cx1.real : dD1 = dD1.copy_negate()
    else : ({}['dD1 Not Recognized.'])
    if cx2.imag == cx1.imag : pass
    elif cx2.imag == -cx1.imag : dD2 = dD2.copy_negate()
    else : ({}['dD2 Not Recognized.'])

    cx3 = complex(dD1, dD2)
    if cx3 == cx1 : return dD1,dD2
    ({}['No match for cx3.'])

make_complex(x)

def make_complex(x, flag = 0) :
    '''
output = make_complex(x [, flag])
Input can be tuple with 1,2 or 3 members.
If 1 or 2 members, 'rect' is understood.
The one member or single object may be int, float, complex, CompleX
or string convertible to int, float or complex.

x = makeComplex(4)
x = makeComplex((4,))
x = makeComplex(('4',0))
x = makeComplex((4,'0', 'rect'))
x = makeComplex(4+0j)
x = makeComplex('4+0j')
x = makeComplex(('4+0j',))
In all seven cases above x = ( Decimal('4'), Decimal('0'), 'rect' )

output is always
(modulus, phase, "polar") or
(real_part, imag_part, "rect")
modulus, phase, real_part, imag_part are Decimal objects.

On error returns None.
'''

    thisName = 'make_complex (x) :'

    if flag :
        print ()
        print (thisName)
        print ('  Input x =', type(x), (str(x))[:60])

    if isinstance(x, CompleX) :  # New class CompleX (note the punctuation.)
        x.check()
        return (x.r, x.i, 'rect')
    
    if isinstance (x,complex) : return make_complex (( x.real, x.imag )) 

    try :
        status = 1
        a,b = str_to_complex (x) 
    except : status = 0
    if status : return a,b,'rect'

    try :
        status = 1
        result = makeDecimal(x)
        isinstance(result, Decimal) or ({}['Expecting result to be type Decimal.'])
        output = result,dD(0),'rect'
    except : status = 0
    if status : return output

    try :
        status = 1
        v1, = x # Allow for (( 3+4j ),)
    except : status = 0
    if status : return make_complex (v1)

    try :
        status = 1 ; error = ''
        if len(x) == 2 :
            v1,v2 = [ makeDecimal (v) for v in x  ]
            str1 = 'rect'
        elif len(x) == 3 :
            v1,v2,str1 = x
            v1,v2 = [ makeDecimal (v) for v in (v1,v2)  ]
        else : ({}['len(x) not in (2,3)'])
        output = v1,v2,str1
        check_complex(output) or ({}['output not valid complex.'])
        if str1[0] in 'rR' : output = v1,v2,'rect'
        else : output = v1,v2,'polar'
    except : status = 0 ; error = str(sys.exc_info())
    if status : return output
    
    if flag : print_error (error)
    return

convertComplex(x)

def convertComplex(x) :
  '''
If input is rectangular, output is polar and vice-versa
'''
  x = make_complex(x)

  if x[2] == 'polar' :
          modulus, phase = x[0], x[1]
          θinRadians = degreesToRadians (phase)
          cosθ = cos(θinRadians)
          sinθ = sin(θinRadians)
          a,b = (modulus*cosθ, modulus*sinθ)
          output = (a,b,'rect')
  else :
      real,imag = x[0], x[1]
      modulus = (real*real + imag*imag).sqrt()
      phase = decATAN2(imag, real)
      output = (modulus, phase, 'polar')

  output = make_complex(output)
  return output

clean_complex (x)

def clean_complex (x) :
    '''
output = clean_complex (x)
output is input returned as complex tuple 
with values "cleaned".
1e-50 is returned as 0.
12.999999999999999..........9999999999999999999999
is returned as 13.
12.888888888888........8888888888888888888888888888 is left unchanged.
Note the following Decimal operations:
>>> getcontext().prec = 20
>>> Decimal('3.9999999999999999999999999999999')
Decimal('3.9999999999999999999999999999999')
>>> Decimal('3.9999999999999999999999999999999')+0
Decimal('4.0000000000000000000')
>>> (Decimal('3.9999999999999999999999999999999')+0).normalize()
Decimal('4')
>>> 
>>> Decimal(76500)
Decimal('76500')
>>> Decimal(76500).normalize()
Decimal('7.65E+4')
>>> Decimal(76500).normalize()+0
Decimal('76500')
>>> 
Hence the line: ((value+0).normalize()+0)
'''
    x = make_complex(x)
    getcontext().prec -= 3

    almostZero = Decimal ('1e-' + str(getcontext().prec) )
    L1 = [
        ( ( (v, Decimal(0))[int(offset)] + 0 ).normalize() + 0 )
        for v in x[:2]
        for offset in (
            (v > -almostZero) and (v < almostZero),
        )
    ]

    for offset in range (0,2,1) :
        v1 = L1[offset]
        if v1 == 0 : pass
        else :
            t1 = v1.as_tuple()
            if len( t1[1] ) < getcontext().prec  : pass
            else : L1[offset] = x[offset]

    getcontext().prec += 3

    while x[2] == 'polar' :
        if L1[0] == 0 :
            L1[1] = L1[0] ; break
        if L1[0] < 0 :
            L1[0] = (L1[0]).copy_negate()
            L1[1] += 180
        L1[1] %= 360
        if L1[1] <= -180 : L1[1] += 360
        if L1[1] > 180 : L1[1] -= 360
        break

    return tuple(L1)

CompleX_to_complex (x)

def CompleX_to_complex (x) :
    '''
output = CompleX_to_complex (x)
output is input formatted for printing.
If x has no imaginary part, output fits on 1 line.
Complex numbers with parts of high precision are returned as string on 2 lines:
( 123_45678_90234_56784_56789_73456_78902_34567_84567.89734_5e-6 -
  345_67890_23456_78456_78973_45678_90234_56782_34567.67654_3e3j )
'''
    def format_Decimal (Decimal1) :
        # This function adds '_' to input string.
        isinstance(Decimal1,str) or ({}['Decimal1 must be type str'])
        if len(Decimal1)<6 : return Decimal1
        str1 = Decimal1.upper() ; sign = ''
        if str1[0] in '+-' :
            sign = str1[0] ; str1 = str1[1:]

        mantissa,exponent = (str1.split('E') + [ '' ])[:2]
        integer_part,decimal_part = (mantissa.split('.') + [ '' ])[:2]
        list1 = list(integer_part)
        for p in range(len(list1)-5, 0, -5) : list1[p:p] = [ '_' ]
        integer_part = ''.join(list1)

        if decimal_part :
            list1 = list(decimal_part)[::-1]
            for p in range(len(list1)-5, 0, -5) : list1[p:p] = [ '_' ]
            decimal_part = ''.join(list1[::-1])
            mantissa = '.'.join((integer_part,decimal_part))
        else : mantissa = integer_part

        if exponent : number = 'E'.join((mantissa,exponent))
        else : number = mantissa

        output = '{}{}'.format(sign,number)
        dD1,dD2 = [ dD(v) for v in (Decimal1,output) ]
        (dD1 == dD2) or ({}['dD1 != dD2'])

        return output

    try :
        error = ''
        x = make_complex(x)
        if x[2] == 'polar' : x = convertComplex (x)
        v1,v2,str1 = x
        (str1 == 'rect') or ({}['str1 must be "rect".'])
        v1,v2 = clean_complex (x)
        if v2 :
            dD1,dD2 = [ dD(str(float(v))) for v in (v1,v2) ]
            set1 = { t1==t2 for t1,t2 in zip((dD1,dD2), (v1,v2)) }
            if set1 == {True} :
                # Each part of complex number converts to float exactly.
                output = complex(dD1,dD2)
            else :
                strr = str(v1) ; stri = str(v2)
                if stri[0] in '+-' :
                    signi = stri[0] ; stri = stri[1:]
                else : signi = '+'
                max = sorted([ len(v) for v in (strr,stri) ])[-1]
                new_line = '''       
     '''.strip(' ')
                joiner = ''
                if max > 50 : joiner = new_line + ' '
                strr = format_Decimal(strr)
                stri = format_Decimal(stri)
                output = '( {} {}{} {}j )'.format(strr,signi,joiner,stri)
        else : output = v1
    except : error = str(sys.exc_info())
    if error :
        list1 = error.split(',')
        print (' ', list1[0])
        print (' ', ','.join(list1[1:-1]) )
        print (' ', list1[-1])
        return
#    sx = 'output' ; print (sx, '=', eval(sx))
    return output

Arithmetic functions

addComplex(v1,v2)

def add_complex(v1,v2) :
    '''
sum = add_complex(v1,v2)
Calculation is rectangular.
The spelling of sum indicates type list or tuple.
Returns None on error.
'''

    try:
        status = 0
        if isinstance (v1,CompleX) : a1,b1 = v1.r,v1.i
        else :
            a1,b1,str1 = make_complex(v1)
            if str1 == 'polar' :
                a1,b1,str1 = convertComplex(v1)
            (str1 == 'rect') or ({}["str1 must be 'rect'."])
        
        if isinstance (v2,CompleX) : a2,b2 = v2.r,v2.i
        else :
            a2,b2,str2 = make_complex(v2)
            if str2 == 'polar' :
                a2,b2,str2 = convertComplex(v2)
            (str2 == 'rect') or ({}["str2 must be 'rect'."])

        sum = (a1+a2, b1+b2)
    except : status = 1
    if status : return
    return sum

def addComplex(v1,v2 = None) :
  '''
SuM = addComplex(list_of_values)
or
SuM = addComplex(v1,v2)
Calculation is rectangular.
The spelling of SuM indicates type CompleX.
Returns None on error.
'''
  try :
      status = 0
      if v2 == None : list_of_values = v1
      else : list_of_values = (v1,v2)
      (len(list_of_values) > 1) or ({}['len(list_of_values) must be > 1'])
      sum = list_of_values[0]
      for v in list_of_values[1:] :
          sum = add_complex(sum,v) 
      
      SuM = CompleX(sum)
      isinstance (SuM, CompleX) or ({}["SuM must be CompleX."])
  except : status = 1
  if status : return
  return SuM

subtractComplex(v1,v2)

def subtractComplex(v1,v2) :
    '''
DifferencE = subtractComplex(v1,v2)
where DifferencE = v1 - v2
Calculation is rectangular.
The spelling of DifferencE indicates type CompleX.
Returns None on error.
'''

    try:
        status = 0
        if isinstance (v1,CompleX) : a,b = v1.r,v1.i
        else :
            a,b,str1 = make_complex(v1)
            if str1 == 'polar' :
                a,b,str1 = convertComplex(v1)
            (str1 == 'rect') or ({}["str1 must be 'rect'."])
        
        if isinstance (v2,CompleX) : c,d = v2.r,v2.i
        else :
            c,d,str2 = make_complex(v2)
            if str2 == 'polar' :
                c,d,str2 = convertComplex(v2)
            (str2 == 'rect') or ({}["str2 must be 'rect'."])

        difference = ( a-c, b-d, 'rect' )
        DifferencE = CompleX(difference)
        isinstance (DifferencE, CompleX) or ({}["DifferencE must be CompleX."])
    except : status = 1
    if status : return
    return DifferencE

multiplyComplex (v1, v2)

def multiply_complex(v1,v2) :
    '''
product = multiply_complex(v1,v2)
Calculation is rectangular.
The spelling of product indicates type list or tuple.
Returns None on error.
'''

    try:
        status = 0
        if isinstance (v1,CompleX) : a,b = v1.r,v1.i
        else :
            a,b,str1 = make_complex(v1)
            if str1 == 'polar' :
                a,b,str1 = convertComplex(v1)
            (str1 == 'rect') or ({}["str1 must be 'rect'."])
        
        if isinstance (v2,CompleX) : c,d = v2.r,v2.i
        else :
            c,d,str2 = make_complex(v2)
            if str2 == 'polar' :
                c,d,str2 = convertComplex(v2)
            (str2 == 'rect') or ({}["str2 must be 'rect'."])

        product = ( a*c-b*d, b*c+a*d, 'rect' )
    except : status = 1
    if status : return
    return product


def multiplyComplex(v1,v2 = None) :
  '''
ProducT = multiplyComplex(list_of_values)
or
ProducT = multiplyComplex(v1,v2)
Calculation is rectangular.
The spelling of ProducT indicates type CompleX.
Returns None on error.
'''
  
  try :
      status = 0
      if v2 == None : list_of_values = v1
      else : list_of_values = (v1,v2)
      (len(list_of_values) > 1) or ({}['len(list_of_values) must be > 1'])
      product = list_of_values[0]
      for v in list_of_values[1:] :
          product = multiply_complex(product,v) 
      
      ProducT = CompleX(product)
      isinstance (ProducT, CompleX) or ({}["ProducT must be CompleX."])
  except : status = 1
  if status : return
  return ProducT

divideComplex (dividend, divisor)

def divideComplex(v1,v2) :
    '''
QuotienT = divideComplex(dividend, divisor)
Calculation is rectangular.
divisor must be non-zero.
if dividend == 0, output is 0.
The spelling of QuotienT indicates type CompleX.
Returns None on error.
'''

    try:
        error = ''
        c,d,str2 = make_complex(v2)
        if str2 == 'polar' : c,d,str2 = convertComplex(v2)
        (str2 == 'rect') or ({}["str2 must be 'rect'."])
        divider = c**2 + d**2
        divider or ({}['divider must be non-zero.'])

        a,b,str1 = make_complex(v1)
        if str1 == 'polar' : a,b,str1 = convertComplex(v1)
        (str1 == 'rect') or ({}["str1 must be 'rect'."])
            
        if a == b == 0 : quotient = 0
        elif d :
            dividendr, dividendi, str3 = multiply_complex( (a,b),(c,-d) )
            (str3 == 'rect') or ({}["str3 must be 'rect'."])
            quotient = dividendr/divider, dividendi/divider, 'rect'
        else : quotient = a/c, b/c, 'rect'

        QuotienT = CompleX(quotient)
        isinstance (QuotienT, CompleX) or ({}["QuotienT must be CompleX."])
    except : error = str(sys.exc_info())
    if error :
        list1 = error.split(',')
        print ('divideComplex(v1,v2) :')
        print (' ', list1[0])
        print (' ', ','.join(list1[1:-1]) )
        print (' ', list1[-1])
        return
    return QuotienT

Exponential functions

complexSQRT (x)

def complexSQRT (x) :
    '''
RooT = complexSQRT (x)
calculation is 'polar'.
'''
    CX1 = CompleX(x)
    if not CX1.m : return CompleX(0)

    modulus, phase = CX1.m, CX1.p
    if modulus < 0 :
            modulus *= -1 ; phase += 180
    modulus = modulus.sqrt()
    phase = phase/2
    root = (modulus,phase,'polar')
    return CompleX(root)

complexCUBEroot (x)

def complexCUBEroot (x) :
    '''
RooT = complexCUBEroot (x)
'polar' output is useful because the other 2 cube roots are:
(root[0], root[1]+120, 'polar')
(root[0], root[1]+240, 'polar')
'''
    CX1 = CompleX(x)                                                                    
    if not CX1.m : return CompleX(0)

# Calculating the cube root is a polar operation.

    modulus, phase = CX1.m, CX1.p
    if modulus < 0 :
            modulus *= -1 ; phase += 180
    modulus_of_root = modulus ** ( Decimal(1) / Decimal(3) )
    phase_of_root = phase/3
    root = (modulus_of_root, phase_of_root, 'polar')
    RooT = CompleX(root)

    return RooT

Functions for testing

cx_sum_zero ( .... )

def cx_sum_zero (values, tolerance=None) :
    '''
sum = cx_sum_zero (values)
If sum is very close to 0, sum is returned as 0.
    '''
    list1 = list(values)
    for p in range(0, len(list1)) :
        v = list1[p]
        if isinstance(v, CompleX) : continue
        list1[p] = CompleX(v)
    
    sum = addComplex(list1)

    max = sorted([ abs(v.m) for v in list1 ])[-1]
    v1 = abs(sum.m)
    if tolerance == None : tolerance = dD('1e-' + str(dgt.prec - 3))
    if v1 < tolerance*max : return CompleX(0)
    return sum

cx_almost_equal (v1,v2)

def cx_almost_equal (v1,v2) :
    '''
status = cx_almost_equal (v1,v2)
status may be : True, False or None
    '''
    def modulus (v) :
        a,b = v
        return (a**2 + b**2).sqrt()

    try :
        error = ''
        v1,v2 = [ make_complex(v) for v in (v1,v2) ]
        a1,b1,str1 = v1
        if str1 == 'polar' : a1,b1,str1 = convertComplex(v1)
        (str1 == 'rect') or ({}['str1 must be "rect".'])
        a2,b2,str2 = v2
        if str2 == 'polar' : a2,b2,str2 = convertComplex(v2)
        (str2 == 'rect') or ({}['str2 must be "rect".'])
        diff = ( a1-a2,b1-b2 ) ; mod1 = modulus(diff)
        sum = ( a1+a2, b1+b2 ) ; mod2 = modulus(sum)
        tolerance = dD('1e-' + str(dgt.prec-4))
        # 
        #                     v1 + v2
        # mod1 <= tolerance * -------
        #                        2
        # 
        status = (2*mod1 <= tolerance*mod2)
    except : error = str(sys.exc_info())
    if error :
        print ('cx_almost_equal (v1,v2) :')
        list1 = error.split(',')
        print (' ', list1[0])
        print (' ', ','.join(list1[1:-1]))
        print (' ', list1[-1])
        sx = 'dgt.prec' ; print (sx, '=', eval(sx))
        return None
    return status

class CompleX

class CompleX :
    '''
This class has 5 attributes:
self.r : the real coordinate of the complex number
self.i : the imaginary coordinate of the complex number
self.m : the modulus of the complex number
self.p : the phase of the complex number in degrees
self.c : the class expressed as Python type complex
'''


    def __init__(self, value=0):
        self.set(value)
        return

    def check(self) :
#        print ('entering check(self) :')
        set1 = { isinstance(v,dD) for v in (self.r, self.i, self.m, self.p) }
        (set1 == {True}) or ({}['Non Decimal in (self.r, self.i, self.m, self.p).'])

        cx1 = CompleX_to_complex ((self.r, self.i, 'rect'))
        if cx1 != self.c : print ("""
For self = {}
Rect values don't match self.c:
cx1    = {}
self.c = {}""".format(
    self,cx1,self.c
            )
)
        cx2 = CompleX_to_complex ((self.m, self.p, 'polar'))
        status = cx_almost_equal (cx2,self.c)
        # status can be True, False, None.
        if not status :
            print ("""
For self = {}
Polar values don't match self.c:
status = {}
cx2    = {}
self.c = {}""".format(
    self,status,cx2,self.c
            )
)
        return

    def set(self, value=Decimal(0)):
#        print ('entering set(self, ',value,'):')
        t1 = make_complex(value)
        if t1[2] == 'rect' :
            self.r, self.i = t1[:2]
            if self.r or self.i :
                t1 = convertComplex(value)
                self.m, self.p = t1[:2]
            else :
                self.m = self.p = Decimal(0)
        else :
            self.m, self.p = t1[:2]
            if self.m :
                t1 = convertComplex(value)
                self.r, self.i = t1[:2]
            else :
                self.r = self.i = self.p = Decimal(0)

        self.c = CompleX_to_complex ((self.r, self.i, 'rect'))
        self.check()
        return

    def clean(self) :
#        print ('entering clean(self) :')
        self.check()
        self.r, self.i = clean_complex ((self.r, self.i, 'rect'))
        self.m, self.p = clean_complex ((self.m, self.p, 'polar'))
        self.check()
        return

    def _print(self) :
        '''
output a string containing all the info
about self and suitable for printing.'''
        self.check()
        return '''
self = {}
real = {}
imag = {}
modulus = {}
phase   = {}
as type complex: {}
'''.format(
    self,
    self.r,    self.i,
    self.m,    self.p,
    self.c,
)

The above code highlights the power of Python in a new class called CompleX. When a new instance of this class is created, it can be accessed through its attributes and changed through its methods with the simplest of syntax.

Seemingly complicated functions such as complex cube root can be implemented in only a few lines of simple code.

Examples

CX1 = CompleX()
print ('isinstance(CX1, CompleX):', isinstance(CX1, CompleX))
print ('CX1 =', CX1._print())

isinstance(CX1, CompleX): True
CX1 =
self = <__main__.CompleX object at 0x101a463c8>
real = 0
imag = 0
modulus = 0
phase = 0
as type complex: 0j

CX2 = CompleX((5,30,'polar'))
print ('isinstance(CX2, CompleX):', isinstance(CX2, CompleX))
print ('CX2 =', CX2._print())

isinstance(CX2, CompleX): True
CX2 =
self = <__main__.CompleX object at 0x101245240>
real = 4.33012701892219323381861585376468
imag = 2.50000000000000000000000000000000
modulus = 5
phase = 30
as type complex: (4.330127018922194+2.5j)

CX3 = CompleX('-5+12j')
print ('isinstance(CX3, CompleX):', isinstance(CX3, CompleX))
print ('CX3 =', CX3._print())

isinstance(CX3, CompleX): True
CX3 =
self = <__main__.CompleX object at 0x101416f28>
real = -5.0
imag = 12.0
modulus = 13.0
phase = 112.619864948040426172949010876680
as type complex: (-5+12j)

CX2 = CX1 # shallow copy.
CX2 = CompleX(CX1) # deep copy.

CX3.check() # This should not produce error.

CX3.set(-7.5)
print ('CX3 =', CX3._print())

self = <__main__.CompleX object at 0x101415f60>
real = -7.5
imag = 0
modulus = 7.5
phase = 180
as type complex: (-7.5+0j)

str1 = '''( - 56876543123452876789012376536455867787653367890876.543246568e-2
            + 12987675645399876553456789023542846223786320765.430987654321E4j )'''
dgt.prec = 50
CX4 = CompleX( str1 )
print ('CX4 =', CX4._print())

CX4 =
self = <__main__.CompleX object at 0x1004ab4d0>
real = -568765431234528767890123765364558677876533678908.76543246568
imag = 129876756453998765534567890235428462237863207654309.87654321
modulus = 1.2987800183682792321609196070852740742247653457504E+50
phase   = 90.250912105532912335540795192437181552898778187515
as type complex: ( -568_76543_12345_28767_89012_37653_64558_67787_65336_78908.76543_24656_8 +
  1_29876_75645_39987_65534_56789_02354_28462_23786_32076_54309.87654_321j )

real, imag retain precision of 60.

modulus, phase have been calculated with precision of 50.

Solving the cubic equation

**Figure 1: Cubic function and X intercepts.**

Origin at point $(0,0)$ .
Intercepts at points $(3,0),\ (-1,0),\ (-2.5,0)$

The cubic equation is expressed thus: $ax^{3}+bx^{2}+cx+d=0$ where both $a,d$ are non-zero.

We will solve the equation $2x^{3}+x^{2}-16x-15=0$ as an example of the use of Python's Decimal module and the new class CompleX. See Figure 1.

dgt.prec = 20

a,b,c,d = (2,1,-16,-15)

A = 2*b*b*b - 9*a*b*c + 27*a*a*d
C = b*b-3*a*c

p = -3*C
q = -A

The depressed cubic and Vieta's substitution

Let $x={\frac {-(b+t)}{3a}}$ , substitute this value of $x$ in the cubic equation and produce the depressed cubic $t^{3}+pt+q=0.$

Let $t=w-{\frac {p}{3w}}$ , substitute this value of $t$ in the depressed equation and produce the quadratic: $W^{2}+qW+C^{3}=0$ where $W=w^{3}$ .

Therefore, $W={\frac {-q\pm {\sqrt {q^{2}-4C^{3}}}}{2}}$ .

disc = q*q - 4*C*C*C

RooT = complexSQRT (disc)

DividenD = addComplex ( -q, RooT )

W = divideComplex(DividenD, 2)

isinstance(W, CompleX) or exit(48)

r3 = Decimal(3).sqrt()

W.clean()

print ('\n################\nW =', W._print())

################
W = 
self = <__main__.CompleX object at 0x104927290>
real = -665
imag = 685.89211979727540825
modulus = 955.33920677422215802
phase   = 134.11396709578581398
as type complex: ( -665 + 685.89211_97972_75408_25j )

The cube roots

**Figure 2: Three cube roots of W.**
$w_{1}$ = RooT1; $w_{2}$ = RooT2; $w_{3}$ = RooT3.
$W\ =\ w_{1}^{3}\ =\ w_{2}^{3}\ =\ w_{3}^{3}$ .
Phase of $W\ =\ \alpha$
Phase of $w_{1}\ =\ \beta \ =\ {\frac {\alpha }{3}}$
Phase of $w_{2}\ =\ \beta \ +\ 120.$
Phase of $w_{3}\ =\ \beta \ -\ 120.$

RooT1 =  complexCUBEroot(W)
t1 = (RooT1.m, RooT1.p+120, 'polar')
t2 = (RooT1.m, RooT1.p-120, 'polar')
RooT2 = CompleX(t1)
RooT3 = CompleX(t2)

for RooT in ( 'RooT1', 'RooT2', 'RooT3' ) :
        root = eval(RooT)
        root.clean()
        print ('++++++++++++++++\n{}:'.format(RooT), root._print())

        print ('''
self.i
------- = {}
sqrt(3)
'''.format(
    float(root.i/r3)
          )
)
# Each value of w has format (k + 1j*m*sqrt(3)). It can be shown that t = 2k.
        t = 2*root.r
        x = -(b+t)/(3*a)
        result = a*x*x*x + b*x*x + c*x + d
        print ('x =', x, 'result =',result)

++++++++++++++++
RooT1:
self = <__main__.CompleX object at 0x104927950>
real = 7
imag = 6.9282032302755091742
modulus = 9.8488578017961047216
phase   = 44.704655698595271327
as type complex: ( 7 + 6.92820_32302_75509_1742j )

self.i
------- = 4.0
sqrt(3)

x = -2.5 result = 0.000
++++++++++++++++
RooT2:
self = <__main__.CompleX object at 0x104927a50>
real = -9.5
imag = 2.5980762113533159395
modulus = 9.8488578017961047216
phase   = 164.70465569859527133
as type complex: ( -9.5 + 2.59807_62113_53315_9395j )

self.i
------- = 1.5
sqrt(3)

x = 3.0 result = 0.000
++++++++++++++++
RooT3:
self = <__main__.CompleX object at 0x104927550>
real = 2.5
imag = -9.5262794416288251140
modulus = 9.8488578017961047216
phase   = -75.295344301404728673
as type complex: ( 2.5 - 9.52627_94416_28825_1140j )

self.i
------- = -5.5
sqrt(3)

x = -1.0 result = 0.000

The three roots of the given cubic are $-1,3,-2.5$ .

Assignments

Calculate the value of $\pi$ to 10,000 places of decimals.

 When calculating  $\pi$ :

>>> Decimal(str(math.pi))
Decimal('3.141592653589793') # π accurate to precision 16.
>>> 
>>> Decimal(math.pi)
Decimal('3.141592653589793115997963468544185161590576171875') # Not an accurate value of π. Why?
>>> #    3.14159265358979323846264338327950288419716939937510582097494459230781 # Correct value of π.

$N$ = $(9,60,$ 'polar' $)$ .

 What are the square roots of  $N$ ?

  $({\sqrt {9}},\ {\frac {60}{2}},$  'polar' $)\ =\ (\pm 3,\ 30,$  'polar' $)$ .  $(-3,\ 30,$  'polar' $)$  =  $(3,\ 30\pm 180,$  'polar' $)$

$r_{1}=(2,60,$ 'polar' $)$ is one of the cube roots of a certain number $N$ .

 What is the number  $N$ ?
 What are the other two cube roots of  $N$ ?

  $r_{1}=2(\cos 60+1j*\sin 60).$ 
 $N=r_{1}^{3}=2^{3}(\cos(60*3)+1j*\sin(60*3))=8(\cos 180+1j*\sin 180)=8(-1)=-8$ .
 The other two cube roots of  $N$  are:  $r_{2},r_{3}=(2,60\pm 120,$  'polar' $)$

Given that CX1 is type CompleX, add method sqrt() to class CompleX so that the expression CX1.sqrt() returns the square root of CX1 as type CompleX.

References

Python's built-in functions:

"abs()",

Python's documentation:

Decimal fixed point and floating point arithmetic, A first look at classes