Re: [boost] [mp_int] new release

Subject: Re: [boost] [mp_int] new release
From: Brandon Kohn (blkohn_at_[hidden])
Date: 2008-10-08 11:45:43

-------------------------------------------------
From: "Kevin Sopp" <baraclese_at_[hidden]>
Sent: Monday, October 06, 2008 8:40 AM
To: <boost_at_[hidden]>
Subject: [boost] [mp_int] new release

> Hi all,
>
> I made a new release of my multiprecision integer library with the
> following changes:
>
> - improved documentation
> - serialization support
> - bugfixes
> - improved benchmark tool
> - prime generation has been reworked
> - probably some other things which I forgot
>
> I believe it's pretty much ready for review now. The mp_int reference
> section is slightly lacking in a few places, but other than that I
> think it's pretty solid.
>

Hi Kevin,

I managed to scrabble together a dumb numeric_limits specialization and was 
able to proceed with my tests. My I'm able to use my lazy_exact type with 
arbitrary precision integer classes by way of boost::rational. So I'm able 
to do comparison tests between your class and mpz_class. There seem to be 
some issues with the mpz_class with respect to the normalization test in 
boost::rational, but I was able to get success by commenting out these 
BOOST_ASSERTs. (not of import to your lib, but wanted to point it out if 
anyone else has tried...)

Anyway, I think I've run into a few bugs with mp_int<> when used in a 
boost::rational.

First there is no A % B operation which is const (boost::rational complained 
as it tries to do a mod operation on a const int_type.) I added one and got 
past this.

During the normalization step of boost::rational the numerator and 
denominator are normalized by the gcd of the two. In one case I was testing 
the %= operator of mp_int<> was causing an infinite loop in the 
boost::math::gcd call. I made an attempt at a fix (though I didn't really 
spend much time on it.. so I'm not sure what if any side-effects)

// TODO currently the sign of the remainder is the sign of the divisor
// i.e. rhs. This is inconsistent with C99, for which it is the sign of
// the dividend. I don't know about C++0x yet.
template<class A, class T>
mp_int<A,T>& mp_int<A,T>::operator %= (const mp_int<A,T>& rhs)
{
  mp_int<A,T> remainder;
  divide(rhs, &remainder);
  //if (remainder.sign_ != rhs.sign_) //! I commented these out to get past 
my issue (as it seems the negative remainder is what the result should be in 
this case)
  //  remainder += rhs;                   //!  I'm sure more scrutiny of the 
sign differences between *this*, rhs, and remainder are required to get it 
right in all cases.
  swap(remainder);
  return *this;
}

The case I had fail was when normalizing (-16384)/(16384). (-1/1). So it was 
infinite looping on boost::math::gcd( mp_int<>( -16384 ), mp_int<>( 
16384 ) );

The next issue I've had which I haven't been able to get past is when I 
initialize an mp_int<> with std::numeric_limits< long >::max(); In this case 
the issue I'm having may not be a bug at all.. I'm not certain. I am 
attempting to create a quick and dirty heuristic to convert a rational into 
a double (for rough calculations in test) and so first tried an iterative 
system which takes an mp_int<> and iteratively divides it by mp_int<>( 
max_long ) until the value will fit into a long. Coupled with some creative 
bookkeeping, my hope is to be able to check if the calculations are in the 
right direction.

Here is what I've done so far:

double rational_cast( const boost::rational< boost::mp_math::mp_int<> >& 
src )
    {
        typedef boost::mp_math::mp_int<> int_t;
        int_t num = src.numerator();
        int_t den = src.denominator();
        int_t maxLong( (std::numeric_limits<long>::max)() );

        int_t q = num / den;
        int_t r = num % den;

        double x=0,n=0,d=0;
        while( q > maxLong )
        {
            x += boost::numeric_cast<double>( 
(std::numeric_limits<long>::max)() );
            q -= maxLong;
        }

        x += boost::numeric_cast<double>( q.to_integral<long>() );

        while( den > maxLong )
        {
            d += boost::numeric_cast<double>( 
(std::numeric_limits<long>::max)() );
            den -= maxLong;
        }

        d += boost::numeric_cast<double>( den.to_integral<long>() );

        while( r > maxLong )
        {
            n += boost::numeric_cast<double>( 
(std::numeric_limits<long>::max)() );
            r -= maxLong;
        }

        n += boost::numeric_cast<double>( r.to_integral<long>() );

        return x + n/d;
    }

I know that there is a possibility that these loops may simply take a very 
large amount of time (and there are probably other issues with what I'm 
trying...). My first calculation with the lazy_exact type using this is to 
approximate pi, and in this case I have confirmed that q is 3. The 
subsequent loops seem to run infinitely however, and inspection of the 
subtractions of maxLong inside mp_int<> do not seem to be fruitful (the 
used_ field never seems to go down, and the value in digits_ bounces around 
significantly large values which do not decrease monotonically as one would 
expect). When I inspected the value for maxLong, digits_ contained the 
correct value, but used_ was set to 1 (is this correct?). The values in n, 
r, and d seemed to have a correspondence between the number of digits 
specified in digits_ field and the value in the used_ field. (I.e. digits_ 
2334974798 has used_ = 10, but digits_ = 2147483647 used_ = 1 for maxLong ).

> I have a few things in mind for the future of this library:
>
> - improved division and squaring algorithms according to the papers on my 
> disk.
> - an implementation of an mp_int proxy type that acts on preallocated
> memory, but does not manage memory itself. This can then be used to
> improve the runtime of several algorithms.
> - expression template support
>
> There is pregenerated html documentation under libs/mp_math/doc/html/
>
> http://www.boostpro.com/vault/index.php?action=downloadfile&filename=mp_math_v03.zip&directory=Math%20-%20Numerics&
>
>
> Kevin Sopp

Respectfully,

Brandon

> _______________________________________________
> Unsubscribe & other changes: 
> http://lists.boost.org/mailman/listinfo.cgi/boost
>