Protoman wrote:
OK, but how do I actually *code* it? The precision alogorithm, I mean.
Before you do the algorithms, settle on a representation. I suggested to use
a std::vector< unsigned long >. For the sake of this discussion, let us go
with that for a while. The idea is that the vector represents a large
number in base 2^32 (I assume for the sake of this discussion that unsigned
long has 32 bits, which may or may not be true for your platform).
Now, you need to implement basic arithmetic for that:
addition,
subtraction,
multiplication,
division
Addition is kind of easy, just start with the least significant digits and
add them. The result is the least significant digit of the sum. If this is
smaller than one of the input digits, a carry occurred. Now move one to the
next digit. Iterate.
Subtraction is essentially the same.
For multiplication, there are choices. You can probably do a straight
forward implementation of the multiplication algorithm that you learned in
elementary school. But that will be slow. Faster methods involve smart
tricks and serious mathematics (including modular arithmetic and discrete
Fourier transform).
Division can be reduced to multiplication. This is tricky but explained in
TAOCP.
Something not to be forgotten: radix conversion in case you want to output
your results. Here is a piece of code that tries to minimize the use of
divisions. It performs reasonably well on something like 1000 digits, but
is very slow compared to libraries available:
namespace DO_NOT_USE {
template < typename T, typename S >
void radix ( T const & num,
std::vector< S > & stack,
std::vector< T > const & power,
typename std::vector< T >::size_type const & power_index,
bool do_fill )
{
if ( power_index == 0 ) {
// one digit only:
stack.push_back ( S( num ) );
} else {
typename std::vector< S >::size_type start = stack.size();
typename std::vector< T >::size_type index = power_index - 1;
T q = num / power[ index ];
T r = num - ( q * power[ index ] );
if ( q == T(0) ) {
radix( r, stack, power, index, false );
} else {
radix( r, stack, power, index, true );
radix( q, stack, power, index, false );
}
if ( do_fill ) {
while ( static_cast<typ ename std::vector< S >::size_type>
( stack.size() - start )
<
static_cast<typ ename std::vector< S >::size_type>
( 1 << power_index ) ) {
stack.push_back ( S(0) );
}
}
}
}
} // namespace DO_NOT_USE
template < typename T, typename S >
std::vector< S > radix ( T const & num, S const & base ) {
std::vector< T > power;
std::vector< S > result;
T current_power = T(base);
T q = num / current_power;
power.push_back ( current_power );
while ( current_power < q ) {
q /= current_power;
current_power *= current_power;
power.push_back ( current_power );
}
T r = num - ( q*current_power );
if ( q == T(0) ) {
DO_NOT_USE::rad ix( r, result, power, power.size()-1, false );
} else {
DO_NOT_USE::rad ix( r, result, power, power.size()-1, true );
DO_NOT_USE::rad ix( q, result, power, power.size()-1, false );
}
return( result );
}
Here T is supposed to be a high-precision cardinal type and S is some other
arithmetic type (like unsigned short). We assume that T can be constructed
from S and converted to S for values in the range of S.
Finally, floating point arithmetic can be reduced to arithmetic of
cardinals.
Also, let me re-iterate that you should use a library. This kind of code is
very time consuming to write, hard to get right, difficult to test, and it
is close to impossible to outperform state of the art implementations .
Best
Kai-Uwe Bux
