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<typename std::vector< S >::size_type>

( stack.size() - start )

<

static_cast<typename 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::radix( r, result, power, power.size()-1, false );

} else {

DO_NOT_USE::radix( r, result, power, power.size()-1, true );

DO_NOT_USE::radix( 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