Re: Large numbers (or is Ada the choice for me?)

["Robert I. Eachus" <rieachus@earthlink.net>]

Hans Georg Schaathun wrote:

> I need a tool to solve large systems of linear equations, with no
> floating point operations (or any other approximations) allowed.
> Even though I am not a seasoned programmer, I think I'll have to
> write the tool myself.
> 
> My question is, will it be reasonably simple to handle large
> rational numbers with Ada?  Is there any packages for this?

Others have answered this correctly, but let me suggest a better way to
approach your original problem.  If the linear equations you are dealing
with are inequalities, then you are trying to solve an (HP-hard) integer
programming problem.  If all of the equations are equalities, then there
is a much easier approach than using abusing a linear algebra package by
instantiating it with a rational type.

Firt find the solution using 64-bit float.  Then, since you stated that
you are only interested in integer solutions, convert the float answers
to integers, and check that all the equations are satisfied.  You may be
able to do this with Ada integer types, but it shouldn't strain the
resouces of the rational arithmetic packages that come with Ada
compilers.

The only remaining problem is that you have to be able to determine
whether or not you are rejecting an otherwise valid solution because the
rounding from float to integer was inaccurate.  If you compute the
determinant of the matrix M in aM = b, and compare it to 0.5/Epsilon you
should be okay.

If your system of equations is too large, or the constants are large,
you may have to resort to the rational arithmetic approach.  An
intermediate choice is to compute the inverse of M above using LUP
decomposition of M.  M' = P'U'L', so since  aM = b, a = bP'U'L'.  P, and
P inverse are permutations of the identity matrix, so they don't add any
error, and U' and L' should be much more accurate than M' itself.  In
fact, by choosing P correctly you should minimize the determinants of U
and L... 
  
> I guess I will manage to implement the rational numbers without
> too much hardship, but I really don't feel like implementing
> arithmetics on large integers.

If you know the size of the maximum intermediates in advance, you can
use the Chinese Remainder trick.  AS long as you choose a vector of
large primes, you can use the extended GCD algorithm to find inverses. 
But always verify that your answers do solve the original equations.  If
any of your intermediate values would be out of range, you will still
get an answer, but possibly an invalid one.