# Solving simultaneous equations

The word "simultaneous" implies that two or more equations exist in the same context at the same time.

For example:

${\displaystyle 3x+2y-13=0\ \dots \ (1)}$

${\displaystyle 2x-y-4=0\ \dots \ (2)}$

Equations ${\displaystyle (1)}$ and ${\displaystyle (2)}$ are 2 simultaneous equations.

When we see two simultaneous equations, it's natural to ask the question: What are the values of ${\displaystyle x}$ and ${\displaystyle y}$ that satisfy both equations ${\displaystyle (1)}$ and ${\displaystyle (2)}$?

To calculate ${\displaystyle (x,y)}$ common to both ${\displaystyle (1)}$ and ${\displaystyle (2)}$:

${\displaystyle (2)*2\ \dots \ 4x-2y-8=0\ \dots \ (3)}$

${\displaystyle (3)+(1)\ \dots \ 7x-21=0;\ x=3.}$

Substitute ${\displaystyle 3}$ for ${\displaystyle x}$ in ${\displaystyle (2):\ 2(3)-y-4=0;\ y=2.}$

The values ${\displaystyle (x,y)=(3,2)}$ satisfy both equations ${\displaystyle (1)}$ and ${\displaystyle (2).}$

In cartesian coordinate geometry of 2 dimensions, point ${\displaystyle (3,2)}$ is the point of intersection of two lines ${\displaystyle (1)}$ and ${\displaystyle (2).}$ See Figure 1 below.

If the 2 lines are parallel, there is no point of intersection. See Figure 2 below.

The trivial solution ${\displaystyle (0,0)}$ is a valid solution. See Figure 3 below.

 Figure 1. 2 lines that intersect at point (3,2). Figure 2. When 2 lines are parallel, there is no point of intersection. Figure 3. The trivial solution x=y=0 is a valid solution.

The information contained in equations ${\displaystyle (1)}$ and ${\displaystyle (2)}$ may constitute a matrix:

```[
[3,  2, -13],
[2, -1,  -4]
]
```

In fact, in the python programming language, the following statement is valid code:

```data = [
[3,  2, -13],
[2, -1,  -4]
]
```

For there to be a point of intersection, `data` must be valid.

Here are examples of invalid data:

```data = [
[3,  2, -13],  # Both lines are parallel.
[6,  4,  -5]   # See Figure 2 above.
]
```
```data = [
[3,  2, -13],
[0,  0,  -4] # Empty row, representing equation 0x + 0y - 4 = 0.
]
```
```data = [
[3,  0, -13], # Empty column representing 2 parallel lines:
[2,  0,  -4]  # 3x = 13 and 2x = 4.
]
```

The matrix above with name "data" is described as matrix "2 by 3," meaning that it contains 2 rows with 3 columns or 3 members per row.

## Solving 2 by 3

Below is python code containing function `solve2by3 (input).` The function accepts as input a matrix like `data` above.

After the code are examples of invocation with both good data and bad data.

```import sys

def solve2by3 (input) :
'''
input represents system:
a1*x + b1*y + c1 = 0
a2*x + b2*y + c2 = 0

input = (
(a1,b1,c1),
(a2,b2,c2),
)

To invoke:
x,y = solve2by3 (input)
Function solve2by3 () may return None.
'''
def solve2by3_(input) :
row1,row2 = input
a1,b1,c1 = row1
a2,b2,c2 = row2

if ( (a1 == b1 == 0) or (a2 == b2 == 0) ) :
print ('solve2by3() : empty row in input.')
return None
if (a1 == a2 == 0) or (b1 == b2 == 0) :
print ('solve2by3() : empty column in input.')
return None
if a1*b2 == a2*b1 :
print ('solve2by3() : parallel input.')
return None
# Arrange input, if necessary, to ensure that a1 is non-zero.
if a1 : pass
else :
row1,row2=row2,row1
a1,b1,c1 = row1
a2,b2,c2 = row2

# The matrix is:
# a1, b1, c1 with a1 non-zero.
# a2, b2, c2
# Process row2, if necessary,  so that a2 becomes 0.

if a2 :
L1 = [a2*v for v in row1]
L2 = [a1*v for v in row2]
zero,b2_,c2_ = [ L1[p]-L2[p] for p in (0,1,2) ]
else : zero,b2_,c2_ =  a2,b2,c2

# The matrix is:
# a1, b1,  c1  with a1 non-zero
# 0,  b2_, c2_ and a2 zero.
# 0*x + b2_*y + c2_ = 0

y = -c2_/b2_
x = -(c1 + b1*y)/a1
return x,y
# Following is python code used for trapping errors
# that may occur in function solve2by3_() above.
output = error = ''
try : output = solve2by3_(input)
except : error = sys.exc_info()[:2]
if (output == None) or error :
if error :
print ('solve2by3() :')
print ('  error detected in processing.')
print ('   ', error)
print ('  input =')
print ('   ', input[0])
print ('   ', input[1])
return None
return output
```

Invocation with good data:

 ```output = solve2by3([ [3, 2, -13], [2, -1, -4] ]) print (output) ``` ```(3.0, 2.0) ```

Invocation with bad or invalid data:

 ```output = solve2by3([ [3, 2, -13], # Both lines are parallel. [6, 4, -5] # See Figure 2 above. ]) print (output) ``` ```solve2by3() : parallel input. input = [3, 2, -13] [6, 4, -5] None ```
 ```output = solve2by3([ [3, 2, -13, 7], # Input matrix not exactly 2 by 3. [6, 4, -5] ]) print (output) ``` ```solve2by3() : error detected in processing. (, ValueError('too many values to unpack (expected 3)',)) input = [3, 2, -13, 7] [6, 4, -5] None ```
 ```output = solve2by3([ [3, 0, 2], [6, 0, 4] ]) print (output) ``` ```solve2by3() : empty column in input. input = [3, 0, 2] [6, 0, 4] None ```

## Solving 3 by 4

Below is python code containing function `solve3by4 (input).` The function accepts as input a 3 by 4 matrix like that described in the code comments.

In 3d coordinate geometry the equations `(1), (2), (3)` are the equations of 3 planes and the solution of the system is the point of intersection of the 3 planes.

Given 3 random planes, there are 2 conditions that do not produce a point of intersection:

• If 2 of the planes are parallel. See Figure 1 below.
• If the 3 planes form a "tent" or "awning." See Figure 2, 3 or 4 below.

The trivial solution ${\displaystyle (0,0,0)}$ is valid. See Figure 5 below.

 Figure 1. When 2 planes are parallel, there is no point of intersection. Figure 2. Triangular tube, 3 planes that never intersect. Figure 3. 3 planes that form a "tent," view 1. There is no point of intersection. Figure 4. 3 planes that form a "tent," view 2. There is no point of intersection. Figure 5. 3 planes that intersect at origin (0, 0, 0).
```def solve3by4 (input) :
'''
input represents system:
a1*x + b1*y + c1*z + d1 = 0 ... (1)
a2*x + b2*y + c2*z + d2 = 0 ... (2)
a3*x + b3*y + c3*z + d3 = 0 ... (3)

input = (
(a1,b1,c1,d1),
(a2,b2,c2,d2),
(a3,b3,c3,d3),
)

To invoke:
x,y,z = solve3by4 (input)
Function solve3by4 () may return None.
'''
def checkDirectionNumbers (input) :
row1,row2,row3 = input
a1,b1,c1,d1 = row1
a2,b2,c2,d2 = row2
a3,b3,c3,d3 = row3
# The direction numbers of the lines of intersection.
dn1 = [b3*c2-b2*c3, a2*c3-a3*c2, b2*a3-b3*a2]
if dn1[0] == dn1[1] == dn1[2] == 0 :
print ('checkDirectionNumbers() : 2 rows parallel, 2 and 3.')
return
dn2 = [b3*c1-b1*c3, a1*c3-a3*c1, b1*a3-b3*a1]
if dn2[0] == dn2[1] == dn2[2] == 0 :
print ('checkDirectionNumbers() : 2 rows parallel, 1 and 3.')
return
dn3 = [b1*c2-b2*c1, a2*c1-a1*c2, b2*a1-b1*a2]
if dn3[0] == dn3[1] == dn3[2] == 0 :
print ('checkDirectionNumbers() : 2 rows parallel, 1 and 2.')
return
# Are 2 lines of intersection parallel?
# If so, and no 2 planes are parallel, the 3 planes form a tent.
a2,b2,c2 = dn2
a3,b3,c3 = dn3
dn1 = [b3*c2-b2*c3, a2*c3-a3*c2, b2*a3-b3*a2]
if dn1[0] == dn1[1] == dn1[2] == 0 :
print ('checkDirectionNumbers() : 3 rows form a tent.');
return

def solve3by4_(input) :
row1,row2,row3 = input
a1,b1,c1,d1 = row1
a2,b2,c2,d2 = row2
a3,b3,c3,d3 = row3

if (a1 == b1 == c1 == 0) or (a2 == b2 == c2 == 0) or (a3 == b3 == c3 == 0) :
print ('solve3by4() : empty row in input.')
return None
if (a1 == a2 == a3 == 0) or (b1 == b2 == b3 == 0) or (c1 == c2 == c3 == 0) :
print ('solve3by4(). empty column in input.')
return None

# Sort input if necessary so that a1 is non-zero.
if a1 : pass
elif a2 :
row1,row2 = row2,row1
a1,b1,c1,d1 = row1
a2,b2,c2,d2 = row2
else :
row1,row3 = row3,row1
a1,b1,c1,d1 = row1
a3,b3,c3,d3 = row3

# The matrix is:
# a1, b1, c1, d1 with a1 non-zero.
# a2, b2, c2, d2
# a3, b3, c3, d3

# Process rows so that a2,a3 become zero.

if a2 :
L1 = [a2*v for v in row1]
L2 = [a1*v for v in row2]
zero,b2_,c2_,d2_ = [ L1[p]-L2[p] for p in (0,1,2,3) ]
else : zero,b2_,c2_,d2_ =  a2,b2,c2,d2

if a3 :
L1 = [a3*v for v in row1]
L2 = [a1*v for v in row3]
zero,b3_,c3_,d3_ = [ L1[p]-L2[p] for p in (0,1,2,3) ]
else : zero,b3_,c3_,d3_ =  a3,b3,c3,d3

# The matrix is:
# a1, b1,  c1,  d1 with a1 non-zero.
# 0,  b2_, c2_, d2_
# 0,  b3_, c3_, d3_
data1 = (
(b2_, c2_, d2_),
(b3_, c3_, d3_),
)
data2 = solve2by3(data1)
if data2 == None :
checkDirectionNumbers (input) # If solve2by3() fails, this line shows why.
return None
y,z = data2
x = -(b1*y + c1*z + d1)/a1

return x,y,z

output = error = ''
try : output = solve3by4_(input)
except : error = sys.exc_info()[:2]

if (output == None) or error :
if error :
print ('solve3by4() :')
print ('  error detected in processing.')
print ('   ', error)
print ('  input of solve3by4() =')
print ('   ', input[0])
print ('   ', input[1])
print ('   ', input[2])
return None
return output
```

Below are examples of invocation with both good data and bad data.

Invocation with good data:

 ```output = solve3by4([ [2, -1, 3, -17], [1, 3, 2, -25] , [3, 2, -1, -12] , ]) print (output) ``` ```(3.0, 4.0, 5.0) ```
 ```output = solve3by4([ [12, -12, -5, -24], [-6, 15, 8, 12] , [ 0, 10, 3, 0] , ]) print (output) ``` ```(2.0, 0.0, -0.0) ```

Invocation with bad or invalid data:

 ```output = solve3by4([ [20, 16, 2, 17], [-77, -27, -25, -24], [81, 45, 18, 43] ]) print (output) ``` ```solve2by3() : parallel input. input = (-692, 346, -829) (396, -198, 517) checkDirectionNumbers() : 3 rows form a tent. input of solve3by4() = [20, 16, 2, 17] [-77, -27, -25, -24] [81, 45, 18, 43] None ```
 ```output = solve3by4([ [-34, -6, 34, -32], [82, 14, -82, 72], [-59, -45, 59, -88] ]) print (output) ``` ```solve2by3() : empty column in input. input = (-16, 0, -176) (-1176, 0, -1104) checkDirectionNumbers() : 3 rows form a tent. input of solve3by4() = [-34, -6, 34, -32] [82, 14, -82, 72] [-59, -45, 59, -88] None ```
 ```output = solve3by4([ [34, -5, -21, 30], [-54, -33, 63, 62], [18, 11, -21, 38] ]) print (output) ``` ```solve2by3() : parallel input. input = (1392, -1008, -3728) (-464, 336, -752) checkDirectionNumbers() : 2 rows parallel, 2 and 3. input of solve3by4() = [34, -5, -21, 30] [-54, -33, 63, 62] [18, 11, -21, 38] None ```
 ```output = solve3by4([ [20, 16, 2, 17], [-77, -27, '-25', -24], [81, 45, 18, 43] ]) print (output) ``` ```solve3by4() : error detected in processing. (, TypeError("unsupported operand type(s) for -: 'int' and 'str'",)) input of solve3by4() = [20, 16, 2, 17] [-77, -27, '-25', -24] [81, 45, 18, 43] None ```

## Solving M by (M+1)

Below is python code containing function `solveMbyN (input).` The function accepts as input a 4 by 5 matrix (or greater) like that described in the code comments.

This function is recursive, calling itself with M being reduced by 1 with each recursive invocation until the matrix is size 3 by 4 in which case function solve3by4() is called.

```def reduceRow (row) :
'''
This function calculates the max(abs) of all values in row,
divides all values by the max(abs) and returns the result.
[27, -15, 42] becomes [27/42, -15/42, 42/42].
Very small values are converted to zero.

To invoke:
output = reduceRow (row)
output is list or None
'''
if False in [ isinstance(v,Decimal) for v in row ] :
print ('reduceRow: non Decimal in input.', { type(v) for v in row })
return None
max = sorted([ abs(v) for v in row ])[-1]
row = [ v/max for v in row ]
almostZero = Decimal('1e-' + str( getcontext().prec ))
row = [ (v, Decimal(0))[abs(v) < almostZero] for v in row ]
return row

def solveMbyN (input) :
'''
input represents system:
a0*w + b0*x + c0*y + d0*z + e0 = 0
a1*w + b1*x + c1*y + d1*z + e1 = 0
a2*w + b2*x + c2*y + d2*z + e2 = 0
a3*w + b3*x + c3*y + d3*z + e3 = 0
or greater.
This solves a matrix of M rows by N columns where N = M+1 and M >= 2

input = [
(a0,b0,c0,d0,e0),
(a1,b1,c1,d1,e1),
(a2,b2,c2,d2,e2),
(a3,b3,c3,d3,e3),
]

To invoke:
w,x,y,z = solveMbyN (input)
Function solveMbyN () may return None.

If input is size 15 by 16:
v1,v2,v3, .... ,v14,v15 = solveMbyN (input)
'''
def solveMbyN_(input) :
numberOfRows = len(input)
numberOfColumns = numberOfRows+1

if numberOfRows < 2 :
print ('solveMbyN()! len(input) < 2.')
return None
if numberOfRows == 2 : return solve2by3(input)
if numberOfRows == 3 : return solve3by4(input)

for row in range (0, numberOfRows) :
if len(input[row]) != numberOfColumns :
print ('solveMbyN()! row with bad length.')
print ('    row number',row,'of',numberOfRows,'rows has length',len(input[row]),'.')
return None
if True not in [ bool(v) for v in input[row][:-1] ] :
print ('solveMbyN()! empty row in input, row number',row,'of',numberOfRows,'rows.')
return None

for column in range (0, numberOfRows) :
sum = 0
for row in range (0, numberOfRows) :
sum = bool(input[row][column])
if sum : break
if not sum :
print ('solveMbyN()! empty column in input. column number',column,'of',numberOfRows,'columns.')
return None

input = list(input)

if 0 : # For matrices greater than 10 by 11 enable this piece of code.
for row in range (0, numberOfRows) :
data = reduceRow (input[row])
if isinstance(data,list) :
input[row] = data

# Arrange input if necessary so that a0 is non-zero.
if not input[0][0] :
for p in range (1, numberOfRows) :
if input[p][0] :
input[p], input[0] = input[0], input[p]
break
# The matrix is:
# [a0,b0,c0,d0,e0], with a0 non-zero,
# [a1,b1,c1,d1,e1],
# [a2,b2,c2,d2,e2],
# [a3,b3,c3,d3,e3],

# Process rows so that a1, a2, a3, .... become zero.
a0 = input[0][0]
for p in range (1, numberOfRows) :
ap = input[p][0]
if ap :
L0 = [ input[0][q]*ap for q in range (0, numberOfColumns) ]
Lp = [ input[p][q]*a0 for q in range (0, numberOfColumns) ]
L_ = [ L0[q]-Lp[q] for q in range (0, numberOfColumns) ]
if L_[0] :
# This should not happen.
print ('solveMbyN()! internal error 1.')
return None
input[p] = L_

# The matrix is:
# [a0, b0, c0, d0, e0], with a0 non-zero,
# [ 0,b1_,c1_,d1_,e1_],
# [ 0,b2_,c2_,d2_,e2_],
# [ 0,b3_,c3_,d3_,e3_],
data1 = [ input[p][1:] for p in range (1, numberOfRows) ]
# data1 = [
# [b1_,c1_,d1_,e1_],
# [b2_,c2_,d2_,e2_],
# [b3_,c3_,d3_,e3_] ]

data2 = solveMbyN(data1) # Here is where function solveMbyN() calls itself.
if data2 == None : return None
topRow = input[0][1:]
# b0, c0, d0, e0 topRow
#  x,  y,  z     data2
L2 = [ topRow[p]*data2[p] for p in range (0, len(data2)) ] + topRow[-1:]
# L2 = [ b0*x, c0*y, d0*z, e0 ]
sum = [ sum for sum in (0,) for p in L2 for sum in (sum+p,) ][-1]
# sum = [ b0*x, b0*x+c0*y, b0*x+c0*y+d0*z, b0*x+c0*y+d0*z+e0 ][-1]
# sum =  b0*x + c0*y + d0*z + e0
# a0*w + b0*x + c0*y + d0*z + e0 = 0
# a0*w + sum = 0

w = -sum/a0

return (w,) + data2

output = error = ''
try : output = solveMbyN_(input)
except : error = sys.exc_info()[:2]
if (output == None) or error :
if error :
print ('solveMbyN()!')
print ('  error detected in processing.')
print ('   ', error)
if len(input) >= 4 :
print ('  input of solveMbyN() with',len(input),'rows =')
for row in input : print ('   ', row)
return None
return output
```

Below are examples of invocation with both good data and bad data.

Invocation with good data:

 ```L1 = [ [Decimal('-870442133162'), Decimal('239132906684'), Decimal('397660499756'), Decimal('788825405004'), Decimal('156830886652'), Decimal('-87096140169')], [Decimal('-112130420352'), Decimal('-93605595828'), Decimal('-185262083472'), Decimal('-93369903383'), Decimal('-497973505065'), Decimal('-766710012513')], [Decimal('-142983492168'), Decimal('619713461197'), Decimal('-608752759383'), Decimal('-448477443042'), Decimal('-281161661751'), Decimal('-799463977245')], [Decimal('572854919541'), Decimal('-148673693908'), Decimal('341291756040'), Decimal('362403188059'), Decimal('-202136092949'), Decimal('-429932229578')], [Decimal('157632861972'), Decimal('330790263276'), Decimal('-836474627955'), Decimal('644295395386'), Decimal('-595590967129'), Decimal('-392835045727')] ] output = solveMbyN(L1) print (output) ``` The following results were produced with decimal precision set at 50: ```( Decimal('-0.06331410760212760088313271882866553044211194694576'), Decimal('1.0969089506763419155148060427588009931309074487534'), Decimal('1.0837038580207682607126722016017867670831200359767'), Decimal('-0.42989566124678842791674322699930461340596671792215'), Decimal('-2.0541601084281682582146753421426610779297051535362') ) ```

Invocation with bad or invalid data:

```output = solveMbyN([
[-4, -6, 3, -4, -8, -8],
[1, 0, 2, 5, 1, 7],
[9, -3, -8, 4, 7, -5],
[1, -9, 5, -2, 5, -2],
[1, -3, -6, -3, -6, 6]
])
print (output)
```
```# With each step towards solve2by3() the length of each number in digits can theoretically double.
solve2by3() : parallel input.
input =
(-181440, 181440, -101088)
(90720, -90720, 406944)
checkDirectionNumbers() : 3 rows form a tent.
input of solve3by4() =
[-756, -1176, 0, -1872]
[-324, -744, 240, -936]
[-324, -384, -120, -264]
input of solveMbyN() with 4 rows =
[-6, 11, 16, -4, 20]
[-66, -5, -20, -44, -92]
[-42, 23, -12, 12, -16]
[-18, -21, -16, -32, 16]
input of solveMbyN() with 5 rows =
[-4, -6, 3, -4, -8, -8]
[1, 0, 2, 5, 1, 7]
[9, -3, -8, 4, 7, -5]
[1, -9, 5, -2, 5, -2]
[1, -3, -6, -3, -6, 6]
None
```

Provided that numbers are small, integers, floats and Decimal objects produce the same results:

 ```output = solveMbyN([ [-4.0, -6.0, 3.0, -4.0, -8.0, -8.0], [1.0, 0.0, 2.0, 5.0, 1.0, 7.0], [9.0, -3.0, -8.0, 4.0, 7.0, -5.0], [1.0, -9.0, 5.0, -2.0, 5.0, -2.0], [1.0, -3.0, -6.0, -3.0, -6.0, 6.0] ]) print (output) ``` ```solve2by3() : parallel input. input = (-181440.0, 181440.0, -101088.0) (90720.0, -90720.0, 406944.0) checkDirectionNumbers() : 3 rows form a tent. input of solve3by4() = [-756.0, -1176.0, 0.0, -1872.0] [-324.0, -744.0, 240.0, -936.0] [-324.0, -384.0, -120.0, -264.0] input of solveMbyN() with 4 rows = [-6.0, 11.0, 16.0, -4.0, 20.0] [-66.0, -5.0, -20.0, -44.0, -92.0] [-42.0, 23.0, -12.0, 12.0, -16.0] [-18.0, -21.0, -16.0, -32.0, 16.0] input of solveMbyN() with 5 rows = [-4.0, -6.0, 3.0, -4.0, -8.0, -8.0] [1.0, 0.0, 2.0, 5.0, 1.0, 7.0] [9.0, -3.0, -8.0, 4.0, 7.0, -5.0] [1.0, -9.0, 5.0, -2.0, 5.0, -2.0] [1.0, -3.0, -6.0, -3.0, -6.0, 6.0] None ```
 ```output = solveMbyN([ [Decimal('-4'), Decimal('-6'), Decimal('3'), Decimal('-4'), Decimal('-8'), Decimal('-8')], [Decimal('1'), Decimal('0'), Decimal('2'), Decimal('5'), Decimal('1'), Decimal('7')], [Decimal('9'), Decimal('-3'), Decimal('-8'), Decimal('4'), Decimal('7'), Decimal('-5')], [Decimal('1'), Decimal('-9'), Decimal('5'), Decimal('-2'), Decimal('5'), Decimal('-2')], [Decimal('1'), Decimal('-3'), Decimal('-6'), Decimal('-3'), Decimal('-6'), Decimal('6')] ]) print (output) ``` ```solve2by3() : parallel input. input = (Decimal('-181440'), Decimal('181440'), Decimal('-101088')) (Decimal('90720'), Decimal('-90720'), Decimal('406944')) checkDirectionNumbers() : 3 rows form a tent. input of solve3by4() = [Decimal('-756'), Decimal('-1176'), Decimal('0'), Decimal('-1872')] [Decimal('-324'), Decimal('-744'), Decimal('240'), Decimal('-936')] [Decimal('-324'), Decimal('-384'), Decimal('-120'), Decimal('-264')] input of solveMbyN() with 4 rows = [Decimal('-6'), Decimal('11'), Decimal('16'), Decimal('-4'), Decimal('20')] [Decimal('-66'), Decimal('-5'), Decimal('-20'), Decimal('-44'), Decimal('-92')] [Decimal('-42'), Decimal('23'), Decimal('-12'), Decimal('12'), Decimal('-16')] [Decimal('-18'), Decimal('-21'), Decimal('-16'), Decimal('-32'), Decimal('16')] input of solveMbyN() with 5 rows = [Decimal('-4'), Decimal('-6'), Decimal('3'), Decimal('-4'), Decimal('-8'), Decimal('-8')] [Decimal('1'), Decimal('0'), Decimal('2'), Decimal('5'), Decimal('1'), Decimal('7')] [Decimal('9'), Decimal('-3'), Decimal('-8'), Decimal('4'), Decimal('7'), Decimal('-5')] [Decimal('1'), Decimal('-9'), Decimal('5'), Decimal('-2'), Decimal('5'), Decimal('-2')] [Decimal('1'), Decimal('-3'), Decimal('-6'), Decimal('-3'), Decimal('-6'), Decimal('6')] None ```

## Solving big matrices

A big matrix is one of size greater than 10 by 11.

With each recursive call of solveMbyN() the size of each value in decimal digits can theoretically double. Python's integer math enjoys infinite precision. If your original matrix is of size 32 by 33, and each value in the matrix is an int of 1 decimal digit, function solve2by3() can be called with ints of size 2**30 (1,073,741,824) decimal digits.

Floats can be used instead of ints, but the limited precision of floats may not provide a desirable level of accuracy in the result.

Python's decimal objects are a good compromise and function reduceRow() above is an attempt to keep the size of each member manageable.

Function reduceRow() works well. However, as with all operations involving python's decimal objects, you need to pay careful attention to precision and accuracy of the results.

On my Mac the following python code generates a matrix of size 100 by 101 and calculates the result containing 100 values in about 1.5 seconds, including the time to generate 10,100 random numbers, each containing up to 12 decimal digits with the decimal point in a random position, and including the time to test output over 100 rows to verify that all calculated values of the solution satisfy each row.

In function solveMbyN() above find line:

```if 0 : # For matrices greater than 10 by 11 enable this piece of code.
```

Change it to:

```if 1 : # For matrices greater than 10 by 11 enable this piece of code.
```

Then:

```from random import getrandbits
from decimal import *
import sys
getcontext().prec=50 # Set precision to 50.

numRows=100

def randomNumber():
# 1 invocation in 256 returns 0.
if not (getrandbits(10) & 0xFF) : return Decimal(0)
return [ w
for q in [ getrandbits(10) & 1]
for numberOfDecimalPlaces in [ getrandbits(10) & 0xF ]
for s1 in [ str(getrandbits(50))[-12:] ]
for s2 in [ '0'*12+s1 ]
for posn in [ len(s2)-numberOfDecimalPlaces ]
for w in [ sign*Decimal(s2[:posn]+'.'+s2[posn:]) ]
][0]

data = [[ randomNumber() for p in range(numRows+1) ]
for p in range (numRows)
]

datan = sorted([ v  for row in data for v in row if v < 0 ])
datap = sorted([ v  for row in data for v in row if v >= 0 ])
print ('A small sample of the random numbers in data:')
print ('+ve values:',len(datap), ',', datap[:2],'\n   ',datap[-2:])
print ('-ve values:',len(datan), ',', datan[:2],'\n   ',datan[-2:])

output = solveMbyN(data)

print ('A small sample of the solutions:')
print ('output[0]  =', output[0])
print ('output[19] =', output[19])
print ('output[59] =', output[59])
print ('output[99] =', output[99])
```
```A small sample of the random numbers in data:
+ve values: 5021 , [Decimal('1.16871265E-7'), Decimal('1.73721245E-7')]
[Decimal('998523464710'), Decimal('999768605679')]
-ve values: 5079 , [Decimal('-991079944209'), Decimal('-989786083110')]
[Decimal('-4.10272824E-7'), Decimal('-2.38664537E-7')]

A small sample of the solutions:
output[0]  = 0.43689400018920695253771832432651757885163362928745
output[19] = 0.19867480862268284151047885891372314311902973127033
output[59] = 0.46169681824826309421940885550155275946668017604646
output[99] = -0.020598220802233733529683040085023443421113627790464
```

### Testing results

Data input to function solveMbyN() is an array of 100 rows, each row containing 101 numbers.

```[
[v0_0, v0_1, v0_2, ..... v0_98,v0_99,v0_100],
.......
[v99_0,v99_1,v99_2, ...... v99_98,v99_99,v99_100],]
]
```

Data received from function solveMbyN() is an array of 100 numbers.

```[ s0,s1,s2, ...... s98,s99 ]
```

The first row is tested with the solutions by;

```L2 = [ v0_0*s0, v0_1*s1, v0_2*s2, ..... v0_98*s98,v0_99*s99,v0_100 ]
```

sum is the sum of all items in L2.

sum is appended to L1 and the process is repeated for all 100 rows.

```L1 = []
for row in data :
# a, b, c, d, e, f, g ....... # row
# u, v, w, x, y, z            # output
L2 = [ row[p]*t1[p] for t1 in (output+(1,),) for p in range (0, len(row)) ]
# L2 = [a*u , b*v , c*w , d*x , e*y , f*z , g]
sum = [ sum for sum in [0] for p in L2 for sum in [sum+p] ][-1]
# sum = [a*u, a*u + b*v, a*u + b*v + c*w, a*u + b*v + c*w + d*x, a*u + b*v + c*w + d*x + e*y,
#        a*u + b*v + c*w + d*x + e*y + f*z, a*u + b*v + c*w + d*x + e*y + f*z + g][-1]
# sum = a*u + b*v + c*w + d*x + e*y + f*z + g
# sum is the sum of all items in L2
L1 += [abs(sum)]

L1a = sorted(L1)
print (len(L1a), L1a[0], L1a[-1])
```

Ideally sum should be zero.

In practice it's extremely unlikely that sum will be zero. My results typically gave values close to

```100 4.1065149609E-38 9.6425614873452398368906E-25
```

The worst of the results above was 23 zeroes after the decimal point.

Better results can be achieved by increasing precision.