VOOZH

URL: https://en.wikiversity.org/wiki/Python_Concepts/Using_Python's_Decimal_module

⇱ Python Concepts/Using Python's Decimal module - Wikiversity

Jump to content

From Wikiversity

< Python Concepts

Objective

[edit | edit source]

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

[edit | edit source]

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

[edit | edit source]

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

>>> fromdecimalimport *
>>> 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 fromcomplex 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

[edit | edit source]

Decimal objects work well with the usual arithmetic operators:

>>> fromdecimalimport *
>>> 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

[edit | edit source]

>>> 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

[edit | edit source]

>>> 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

[edit | edit source]

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

[edit | edit source]

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:

>>> fromdecimalimport *
>>>
>>> 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.

>>> fromdecimalimport *
>>>
>>> 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.

>>> fromdecimalimport *
>>> 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():

>>> fromdecimalimport *
>>>
>>> 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:

>>> fromdecimalimport *
>>> 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

[edit | edit source]

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

[edit | edit source]

defprint_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.

importdecimal

getcontext = decimal.getcontext

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

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

defmakeDecimal (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

[edit | edit source]

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

For the functions 👁 {\displaystyle \cos x}
and 👁 {\displaystyle \sin x}
see Recipes.

arctan x

[edit | edit source]

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

👁 {\displaystyle \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}

defarctan (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 π

[edit | edit source]

There are many ways to calculate 👁 {\displaystyle \pi }
. For example:

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

Based on the following, we will perform six calculations of 👁 {\displaystyle \pi }
and compare the results.

👁 {\displaystyle \pi =30*\arctan(\tan 6)}

👁 {\displaystyle \pi =24*\arctan(\tan 7.5)}

👁 {\displaystyle \pi =12*\arctan(\tan 15)}

👁 {\displaystyle \pi =8*\arctan(\tan 22.5)}

👁 {\displaystyle \pi =6*\arctan(\tan 30)}

👁 {\displaystyle \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)

[edit | edit source]

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

deftanθ_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)

[edit | edit source]

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

This function invokes 👁 {\displaystyle arctan(x)\ 0\leq x\leq 1}
.

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

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

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

👁 {\displaystyle 68=4\arctan(\tan {\frac {68}{4}})}

👁 {\displaystyle 32=2\arctan(\tan {\frac {32}{2}})}

👁 {\displaystyle 12=\arctan(\tan 12)}

defdecATAN2 (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)

[edit | edit source]

defdegreesToRadians (θ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

[edit | edit source]

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 👁 {\displaystyle modulus\ *\ (\ \cos(phaseInDegrees)\ +\ 1j*\sin(phaseInDegrees)\ )}
.

Or:

(realPart, imaginaryPart, 'rect')

where 👁 {\displaystyle realPart,imaginaryPart}
are the rectangular coordinates of complex number 👁 {\displaystyle realPart\ +\ 1j*imaginaryPart}
.

The four values 👁 {\displaystyle 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:

👁 {\displaystyle -n*(\cos \phi +1j*\sin \phi )=n*(\cos(\phi \pm 180)+1j*\sin(\phi \pm 180))}

and, for example:

👁 {\displaystyle {\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 👁 {\displaystyle -4}
is:

👁 {\displaystyle {\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

[edit | edit source]

check_complex(x)

[edit | edit source]

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

defcheck_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)

[edit | edit source]

This function accepts a string like 👁 {\displaystyle (2345678901234567.8903456789-9876543219654321098.765432765j)}
and preserves precision of both real and imaginary parts.

importre

# 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)

defstr_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)

[edit | edit source]

defmake_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)

[edit | edit source]

defconvertComplex(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)

[edit | edit source]

defclean_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)

[edit | edit source]

defCompleX_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 )
'''
 defformat_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

[edit | edit source]

addComplex(v1,v2)

[edit | edit source]

defadd_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

defaddComplex(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)

[edit | edit source]

defsubtractComplex(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)

[edit | edit source]

defmultiply_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


defmultiplyComplex(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)

[edit | edit source]

defdivideComplex(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

[edit | edit source]

complexSQRT (x)

[edit | edit source]

defcomplexSQRT (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)

[edit | edit source]

defcomplexCUBEroot (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

[edit | edit source]

cx_sum_zero ( .... )

[edit | edit source]

defcx_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)

[edit | edit source]

defcx_almost_equal (v1,v2) :
'''
status = cx_almost_equal (v1,v2)
status may be : True, False or None
 '''
 defmodulus (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

[edit | edit source]

classCompleX :
'''
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

 defcheck(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

 defset(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

 defclean(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

[edit | edit source]

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

[edit | edit source]

Figure 1: Cubic function and X intercepts.

Origin at point 👁 {\displaystyle (0,0)}
.
Intercepts at points 👁 {\displaystyle (3,0),\ (-1,0),\ (-2.5,0)}

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

We will solve the equation 👁 {\displaystyle 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

[edit | edit source]

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

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

Therefore, 👁 {\displaystyle 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

[edit | edit source]

Figure 2: Three cube roots of W.
👁 {\displaystyle w_{1}}
= RooT1; 👁 {\displaystyle w_{2}}
= RooT2; 👁 {\displaystyle w_{3}}
= RooT3.
👁 {\displaystyle W\ =\ w_{1}^{3}\ =\ w_{2}^{3}\ =\ w_{3}^{3}}
.
Phase of 👁 {\displaystyle W\ =\ \alpha }

Phase of 👁 {\displaystyle w_{1}\ =\ \beta \ =\ {\frac {\alpha }{3}}}

Phase of 👁 {\displaystyle w_{2}\ =\ \beta \ +\ 120.}

Phase of 👁 {\displaystyle 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 👁 {\displaystyle -1,3,-2.5}
.

Assignments

[edit | edit source]

Calculate the value of 👁 {\displaystyle \pi }
to 10,000 places of decimals.

 When calculating 👁 {\displaystyle \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 π.

👁 {\displaystyle N}
= 👁 {\displaystyle (9,60,}
'polar'👁 {\displaystyle )}
.

 What are the square roots of 👁 {\displaystyle N}
?

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

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

 What is the number 👁 {\displaystyle N}
?
 What are the other two cube roots of 👁 {\displaystyle N}
?

 👁 {\displaystyle r_{1}=2(\cos 60+1j*\sin 60).}

👁 {\displaystyle 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 👁 {\displaystyle N}
 are: 👁 {\displaystyle r_{2},r_{3}=(2,60\pm 120,}
 'polar'👁 {\displaystyle )}

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.

Further Reading or Review

[edit | edit source]

References

[edit | edit source]

Python's built-in functions:

Python's documentation:

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

Retrieved from "https://en.wikiversity.org/w/index.php?title=Python_Concepts/Using_Python%27s_Decimal_module&oldid=2745685"

Python

Hidden category:

Pages with syntax highlighting errors