| > using arbitrary precision fractions for Gaussian elimination Here's a way to get exact results from just standard precision, say, integers 32 bits long: For Gauss elimination, that's for solving a system of linear equations, say, m equations in n unknowns. So, we are given, in floating point, a matrix m x n A, an unknown vector 1 x n x, and a right side constant 1 x m b. So we are looking for x so that Ax = b Now, multiply each equation by an integer, say, a power of 2, so that all the numbers in A and b are integers. Next get a list of positive prime numbers, each, say, stored in an integer 32 bits long. Say we have p such prime numbers; so we have prime number for i = 1, 2, ..., p. Next for each i, solve Ax = b by using Gauss elimination but in the field of integers modulo prime number i. For division, use the Euclidean greatest common divisor algorithm. Yes for this arithmetic have to be able to form the 64 bit sum or product of two 32 bit whole numbers and then divide by 32 bit integer and keep the 32 bit remainder -- commonly we write a little assembler routine for this. After doing that work p times (could be done in parallel), we use the Chinese remainder theorem to put together the rational numbers that are quotients of multi-precision whole numbers of the solution. With those, we can get floating point approximations as accurate as we please. But if want to work in floating point arithmetic anyway, then there are three ways to do better: (1) To start, in Gauss elimination, look down column 1 of A, find the row with the number of largest absolute value, and swap rows to put that number in row 1 and column 1. More generally after p - 1 rows, will want a pivot element in row p and column p. By swapping rows, use the number largest in absolute value in column p and rows p, p + 1, ..., m. (2) When form an inner product, say, u(1)v(1) + u(2)v(2) + ... + u(q)v(q) form each product in double precision and add the double precision numbers. If want to do a little better, then before adding, sort the numbers so that are adding the smaller numbers first. Or "The main sin in numerical analysis is subtracting two large numbers whose difference is small.", etc. (3) Once have x so that Ax ~ b, find difference d = b - Ax and then solve for e in Ae = d and replace x by x + e so that A(x + e) = Ax + Ae = (b - d) + d = b. So, take x + e as the solution. After Gauss elimination, solving Ae = d can go relatively quickly. Can do this step more than once. |