sobota 19. prosince 2009

Jde to i přesně

Problém zaokrouhlování je nejen nepřesnost, ale i neplatnost asociativního zákona i když komutativita zůstává neporušena. Např:

(-1010 + 10-10) + 1010 = -1010 + 1010 = 0 ≠ (-1010 + 1010) + 10-10 = 10-10.

Výsledek je tedy nejen nepřesný, ale i pokaždé jiný v závislosti na tom, v jakém pořadí jsou operace vykonávány. Obě nevýhody lze zcela odstranit, pokud se omezíme jen na racionální čísla. To lze v případě řešení soustavy lineárních rovnic, protože řešení soustavy s racionálními koeficienty je opět racionální, což lze snadno dokázat. Přiložený skript je implementace Gaussovy eliminační metody. Řešení (ovšem až na výpis) je zcela přesné bez jakéhokoliv zaokrouhlování, což je zaručeno tím, že racionální číslo je reprezentováno dvěmi čísly celými. Skriptovací jazyk Python nemá kromě omezené paměti jiné omezení na celá čísla.


def codeInt(x):
return [x, 1]

def RatMPrint(A):
for r in A:
for c in r:
print float(c[0]) / float(c[1]),
print

def RatEq(a, b):
return a[0] * b[1] == b[0] * a[1]

def RatAdd(a, b):
return [a[0] * b[1] + b[0] * a[1], a[1] * b[1]]

def RatDiff(a, b):
return [a[0] * b[1] - b[0] * a[1], a[1] * b[1]]

def RatMult(a, b):
return [a[0] * b[0], a[1] * b[1]]

def RatDiv(a, b):
return [a[0] * b[1], a[1] * b[0]]

def RatAbs(a):
return [abs(a[0]), abs(a[1])]

def RatIsGr(a, b):
return a[0] * b[1] > b[0] * a[1]

def RatSize(A):
r = len(A)
c = 0
if r > 0:
c = len(A[0])
return [r, c]

def RatMGaussJordan(m):
[h, w] = RatSize(m)
for y in range(0, h):
maxrow = y
for y2 in range(y + 1, h): # Find max pivot
if RatIsGr(RatAbs(m[y2][y]), RatAbs(m[maxrow][y])):
maxrow = y2
(m[y], m[maxrow]) = (m[maxrow], m[y])
if RatEq(m[y][y], codeInt(0)): # Singular?
return False
for y2 in range(y + 1, h): # Eliminate column y
c = RatDiv(m[y2][y], m[y][y])
for x in range(y, w):
m[y2][x] = RatDiff(m[y2][x], RatMult(m[y][x], c))
for y in range(h-1, 0-1, -1): # Backsubstitute
c = m[y][y]
for y2 in range(0, y):
for x in range(w - 1, y - 1, -1):
m[y2][x] = RatDiff(m[y2][x], RatMult(m[y][x], RatDiv(m[y2][y], c)))
m[y][y] = RatDiv(m[y][y], c)
for x in range(h, w): # Normalize row y
m[y][x] = RatDiv(m[y][x], c)
return True

#*******************************************************************************

m = [[codeInt(1), codeInt(2), codeInt(0)],
[codeInt(3), codeInt(4), codeInt(1)]]
if RatMGaussJordan(m):
RatMPrint(m)