From: martin@math.rwth-aachen.de (  Martin Schoenert)
Newsgroups: comp.programming
Subject: Re: Arbitrary precision with large bases
Date: 6 Jan 93 09:10:24 GMT
Organization: Rechnerbetrieb Informatik - RWTH Aachen
NNTP-Posting-Host: bert.math.rwth-aachen.de

markh@csd4.csd.uwm.edu (Mark) writes:

       A while back, I came up with an efficient arbitrary precision BCD
    arithmetic package that uses base 1,000,000,000 (which happens to be
    close to 2^30).

       I'm looking for other packages for comparison that use large bases,
    either BCD or binary.  The non-trivial aspects to such as package are
    the following:

       (1) An algorithm to carry out A, B -> AB div BASE, AB mod BASE
       (2) A division algorithm that minimizes the amount of backtracking.
           The grade school algorithm is a backtracking algorithm, and only
           in base 2 is backtracking avoided.
       (3) A square root algorithm, like the one taught in school, that likewise
           minimizes backtracking.

I have looked at quite a few large integer arithmetic packages and have
written two myself.  Mark Riordan has collected several, they are
available for anonymous 'ftp' on 'cl-next2.cl.msu.edu' (35.8.4.22).

Most packages use a base that is a power of 2.  The most common choices
are 2^16 and 2^32.  I assume that 2^16 does not count as a large base
because with base 2^16 computing AB div BASE and AB mod BASE is not a
problem at all.

Base 2^32 has the following advantages.

a)  the larger the base the fewer digits any given integer has in this
    base.  Since (for example) the performance of the multiplication
    is proportional to the product of the number of digits of the
    operands it is always a win to use larger bases (though the win
    with base 2^32 over base 10^9 is only about 12.5%).

b)  many modern processors have a 32bit x 32bit -> 64bit instruction.
    This instruction cannot directly be used from within C, so one needs
    a little bit of assembler.  However if one does that one gets
    AB div BASE and AB mod BASE with one instruction.

Of course base 2^32 has the disadvantage that converting an integer to
human readable format (i.e., base 10) is more work than in base 10^9.

With respect to aspect (2).  Knuth shows that if you normalize the
divisor such that its leading digit is larger than or equal to BASE/2 and
compute an appoximation of a digit of the quotient by dividing the
leading three digits of the remainder by the leading two digits of the
divisor the approximation is at most too large by 1 *independent* of the
base.  It is easy to detect if the digit is too large and make the
appropriate corrections.

With respect to aspect (3).  The  usual method to compute square roots is
Newton's.

Here is another method, which I think is just a careful implementation of
what you refer to as the scool method.  This method is considerably
faster than Newton's, in principle it is as fast as squaring the result,
in practice it is about a factor of 2 slower in my arithmetic.

Let BB be a power of BASE and let

    A  =  A_1 BB^2 + A_0

such that BB^2/4 <= A_1.  Compute the upper half of the root recursively

    X_1  =  Floor( Root( A_1 ) );
    R  =  A - (X_1 BB)^2;

now compute the approximation for the rest of the root

    X_0  =  R / (2 X_1 BB);
    S  =  R - 2 X_0 X_1 BB;

this approximation may be off by one, correct if neccessary

    if S < X_0^2  then
        X_0  =  X_0 - 1;
        S  =  S + 2 X_1 BB;  /* = R - 2 X_0 X_1 BB again */
    fi;

then with

    X  =  X_1 BB + X_0;
    T  =  S - X_0^2;

we have

    X  =  Floor( Root( A ) );
    T  =  A - X^2;

You can either choose BB to be BASE, in which case each iteration
computes one more digit of the root, or you can choose BB such that A_1
and A_0 have the same number of digits, in which case each iteration
doubles the number of digits of the root.  The second is faster for very
large integers (~1000 decimal digits) if one uses Karatsuba's trick for
the computation of X_0^2.

This algorithm is again independent of the value of BASE.

Hope this helps,  Martin.

-- .- .-. - .. -.  .-.. --- ...- . ...  .- -. -. .. -.- .-
Martin Sch"onert,   Martin.Schoenert@Math.RWTH-Aachen.DE,  +49 241 804551
Lehrstuhl D f"ur Mathematik, Templergraben 64, RWTH, D 51 Aachen, Germany