source: git/factory/int_rat.cc @ ed66770

/* emacs edit mode for this file is -*- C++ -*- */
/* $Id$ */

#include "config.h"

#include "cf_assert.h"

#include "cf_defs.h"
#include "int_rat.h"
#include "int_int.h"
#include "imm.h"
#include "canonicalform.h"
#include "cf_factory.h"
#include "gmpext.h"

static int intgcd( int a, int b )
{
    if ( a < 0 ) a = -a;
    if ( b < 0 ) b = -b;

    int c;

    while ( b != 0 )
    {
        c = a % b;
        a = b;
        b = c;
    }
    return a;
}


InternalRational::InternalRational()
{
    mpz_init( _num );
    mpz_init_set_si( _den, 1 );
}

InternalRational::InternalRational( const int i )
{
    mpz_init_set_si( _num, i );
    mpz_init_set_si( _den, 1 );
}

InternalRational::InternalRational( const int n, const int d )
{
    ASSERT( d != 0, "divide by zero" );
    if ( n == 0 )
    {
        mpz_init_set_si( _num, 0 );
        mpz_init_set_si( _den, 1 );
    }
    else
    {
        int g = intgcd( n, d );
        if ( d < 0 )
        {
          mpz_init_set_si( _num, -n / g );
          mpz_init_set_si( _den, -d / g );
        }
        else
        {
          mpz_init_set_si( _num, n / g );
          mpz_init_set_si( _den, d / g );
        }
    }
}

InternalRational::InternalRational( const char * )
{
    // sollte nicht gebraucht werden !!!
    ASSERT( 0, "fatal error" );
    mpz_init( _num );
    mpz_init( _den );
}

//InternalRational::InternalRational( const mpz_ptr n ) : _num(n)
//{
//    mpz_init_set_si( _den, 1 );
//}

InternalRational::InternalRational( const mpz_ptr n )
{
    _num[0]=*n;
    mpz_init_set_si( _den, 1 );
}

InternalRational::InternalRational( const mpz_ptr n, const mpz_ptr d )
{
  _num[0]=*n;
  _den[0]=*d;
}

InternalRational::~InternalRational()
{
    mpz_clear( _num );
    mpz_clear( _den );
}

InternalCF* InternalRational::deepCopyObject() const
{
    mpz_t dummy_num;
    mpz_t dummy_den;
    mpz_init_set( dummy_num, _num );
    mpz_init_set( dummy_den, _den );
    return new InternalRational( dummy_num, dummy_den );
}

#ifndef NOSTREAMIO
void InternalRational::print( OSTREAM & os, char * c )
{
    char * str = new char[mpz_sizeinbase( _num, 10 ) + 2];
    str = mpz_get_str( str, 10, _num );
    os << str << '/';
    delete [] str;
    str = new char[mpz_sizeinbase( _den, 10 ) + 2];
    str = mpz_get_str( str, 10, _den );
    os << str << c;
    delete [] str;
}
#endif /* NOSTREAMIO */

bool InternalRational::is_imm() const
{
    return mpz_cmp_si( _den, 1 ) == 0 && mpz_is_imm( _num );
}

InternalCF* InternalRational::genZero()
{
    if ( isZero() )
        return copyObject();
    else
        return new InternalRational();
}

InternalCF* InternalRational::genOne()
{
    if ( isOne() )
        return copyObject();
    else
        return new InternalRational( 1 );
}

//{{{ InternalCF * InternalRational::num (), den ()
// docu: see CanonicalForm::num(), den()
InternalCF * InternalRational::num ()
{
    if ( mpz_is_imm( _num ) )
    {
        InternalCF * res = int2imm( mpz_get_si( _num ) );
        return res;
    }
    else
    {
        mpz_t dummy;
        mpz_init_set( dummy, _num );
        return new InternalInteger( dummy );
    }
}

InternalCF * InternalRational::den ()
{
    if ( mpz_is_imm( _den ) )
    {
        InternalCF * res = int2imm( mpz_get_si( _den ) );
        return res;
    }
    else
    {
        mpz_t dummy;
        mpz_init_set( dummy, _den );
        return new InternalInteger( dummy );
    }
}
//}}}

//{{{ InternalCF * InternalRational::neg ()
// docu: see CanonicalForm::operator -()
InternalCF *
InternalRational::neg ()
{
    if ( getRefCount() > 1 )
    {
        decRefCount();
        mpz_t dummy_num;
        mpz_t dummy_den;
        mpz_init_set( dummy_num, _num );
        mpz_init_set( dummy_den, _den );
        mpz_neg( dummy_num, dummy_num );
        return new InternalRational( dummy_num, dummy_den );
    }
    else
    {
        mpz_neg( _num, _num );
        return this;
    }
}
//}}}

InternalCF* InternalRational::addsame( InternalCF * c )
{
    ASSERT( ! ::is_imm( c ) && c->levelcoeff() == RationalDomain, "illegal domain" );
    mpz_t n, d, g;

    mpz_init( g ); mpz_init( n ); mpz_init( d );
    mpz_gcd( g, _den, MPQDEN( c ) );

    if ( mpz_cmp_si( g, 1 ) == 0 )
    {
      mpz_mul( n, _den, MPQNUM( c ) );
      mpz_mul( g, _num, MPQDEN( c ) );
      mpz_add( n, n, g );
      mpz_mul( d, _den, MPQDEN( c ) );
    }
    else
    {
      mpz_t tmp1;
      mpz_t tmp2;
      mpz_init( tmp1 );
      mpz_divexact( tmp1, _den, g );
      mpz_init( tmp2 );
      mpz_divexact( tmp2, MPQDEN( c ), g );
      mpz_mul( d, tmp2, _den );
      mpz_mul( tmp2, tmp2, _num );
      mpz_mul( tmp1, tmp1, MPQNUM( c ) );
      mpz_add( n, tmp1, tmp2 );
      mpz_gcd( g, n, d );
      if ( mpz_cmp_si( g, 1 ) != 0 )
      {
          mpz_divexact( n, n, g );
          mpz_divexact( d, d, g );
      }
      mpz_clear( tmp1 );
      mpz_clear( tmp2 );
    }
    mpz_clear( g );
    if ( deleteObject() ) delete this;
    if ( mpz_cmp_si( d, 1 ) == 0 )
    {
      mpz_clear( d );
      if ( mpz_is_imm( n ) )
      {
          InternalCF * res = int2imm( mpz_get_si( n ) );
          mpz_clear( n );
          return res;
      }
      else
      {
          return new InternalInteger( n );
      }
    }
    else
    {
      return new InternalRational( n, d );
    }
}

InternalCF* InternalRational::subsame( InternalCF * c )
{
    ASSERT( ! ::is_imm( c ) && c->levelcoeff() == RationalDomain, "illegal domain" );
    mpz_t n, d, g;

    mpz_init( g ); mpz_init( n ); mpz_init( d );
    mpz_gcd( g, _den, MPQDEN( c ) );

    if ( mpz_cmp_si( g, 1 ) == 0 )
    {
        mpz_mul( n, _den, MPQNUM( c ) );
        mpz_mul( g, _num, MPQDEN( c ) );
        mpz_sub( n, g, n );
        mpz_mul( d, _den, MPQDEN( c ) );
    }
    else
    {
        mpz_t tmp1;
        mpz_t tmp2;
        mpz_init( tmp1 );
        mpz_divexact( tmp1, _den, g );
        mpz_init( tmp2 );
        mpz_divexact( tmp2, MPQDEN( c ), g );
        mpz_mul( d, tmp2, _den );
        mpz_mul( tmp2, tmp2, _num );
        mpz_mul( tmp1, tmp1, MPQNUM( c ) );
        mpz_sub( n, tmp2, tmp1 );
        mpz_gcd( g, n, d );
        if ( mpz_cmp_si( g, 1 ) != 0 )
        {
            mpz_divexact( n, n, g );
            mpz_divexact( d, d, g );
        }
        mpz_clear( tmp1 );
        mpz_clear( tmp2 );
    }
    mpz_clear( g );
    if ( deleteObject() ) delete this;
    if ( mpz_cmp_si( d, 1 ) == 0 )
    {
        mpz_clear( d );
        if ( mpz_is_imm( n ) )
        {
            InternalCF * res = int2imm( mpz_get_si( n ) );
            mpz_clear( n );
            return res;
        }
        else
        {
            return new InternalInteger( n );
        }
    }
    else
        return new InternalRational( n, d );
}

InternalCF* InternalRational::mulsame( InternalCF * c )
{
    ASSERT( ! ::is_imm( c ) && c->levelcoeff() == RationalDomain, "illegal domain" );
    mpz_t n, d;
    mpz_init( n ); mpz_init( d );

    if ( this == c )
    {
        mpz_mul( n, _num, _num );
        mpz_mul( d, _den, _den );
    }
    else
    {
        mpz_t g1, g2, tmp1, tmp2;
        mpz_init( g1 ); mpz_init( g2 );
        mpz_gcd( g1, _num, MPQDEN( c ) );
        mpz_gcd( g2, _den, MPQNUM( c ) );
        bool g1is1 = mpz_cmp_si( g1, 1 ) == 0;
        bool g2is1 = mpz_cmp_si( g2, 1 ) == 0;
        mpz_init( tmp1 ); mpz_init( tmp2 );
        if ( ! g1is1 )
            mpz_divexact( tmp1, _num, g1 );
        else
            mpz_set( tmp1, _num );
        if ( ! g2is1 )
            mpz_divexact( tmp2, MPQNUM( c ), g2 );
        else
            mpz_set( tmp2, MPQNUM( c ) );
        mpz_mul( n, tmp1, tmp2 );
        if ( ! g1is1 )
            mpz_divexact( tmp1, MPQDEN( c ), g1 );
        else
            mpz_set( tmp1, MPQDEN( c ) );
        if ( ! g2is1 )
            mpz_divexact( tmp2, _den, g2 );
        else
            mpz_set( tmp2, _den );
        mpz_mul( d, tmp1, tmp2 );
        mpz_clear( tmp1 ); mpz_clear( tmp2 );
        mpz_clear( g1 ); mpz_clear( g2 );
    }
    if ( deleteObject() ) delete this;
    if ( mpz_cmp_si( d, 1 ) == 0 )
    {
        mpz_clear( d );
        if ( mpz_is_imm( n ) )
        {
            InternalCF * res = int2imm( mpz_get_si( n ) );
            mpz_clear( n );
            return res;
        }
        else
        {
            return new InternalInteger( n );
        }
    }
    else
        return new InternalRational( n, d );
}

InternalCF* InternalRational::dividesame( InternalCF * c )
{
    ASSERT( ! ::is_imm( c ) && c->levelcoeff() == RationalDomain, "illegal domain" );

    if ( this == c )
    {
        if ( deleteObject() ) delete this;
        return CFFactory::basic( 1 );
    }
    else
    {
        mpz_t n, d;
        mpz_t g1, g2, tmp1, tmp2;
        mpz_init( n ); mpz_init( d );
        mpz_init( g1 ); mpz_init( g2 );
        mpz_gcd( g1, _num, MPQNUM( c ) );
        mpz_gcd( g2, _den, MPQDEN( c ) );
        bool g1is1 = mpz_cmp_si( g1, 1 ) == 0;
        bool g2is1 = mpz_cmp_si( g2, 1 ) == 0;
        mpz_init( tmp1 ); mpz_init( tmp2 );
        if ( ! g1is1 )
            mpz_divexact( tmp1, _num, g1 );
        else
            mpz_set( tmp1, _num );
        if ( ! g2is1 )
            mpz_divexact( tmp2, MPQDEN( c ), g2 );
        else
            mpz_set( tmp2, MPQDEN( c ) );
        mpz_mul( n, tmp1, tmp2 );
        if ( ! g1is1 )
            mpz_divexact( tmp1, MPQNUM( c ), g1 );
        else
            mpz_set( tmp1, MPQNUM( c ) );
        if ( ! g2is1 )
            mpz_divexact( tmp2, _den, g2 );
        else
            mpz_set( tmp2, _den );
        mpz_mul( d, tmp1, tmp2 );
        mpz_clear( tmp1 ); mpz_clear( tmp2 );
        mpz_clear( g1 ); mpz_clear( g2 );
        if ( deleteObject() ) delete this;
        if ( mpz_cmp_si( d, 0 ) < 0 )
        {
            mpz_neg( d, d );
            mpz_neg( n, n );
        }
        if ( mpz_cmp_si( d, 1 ) == 0 )
        {
            mpz_clear( d );
            if ( mpz_is_imm( n ) )
            {
                InternalCF * res = int2imm( mpz_get_si( n ) );
                mpz_clear( n );
                return res;
            }
            else
            {
                return new InternalInteger( n );
            }
        }
        else
            return new InternalRational( n, d );
    }
}

InternalCF* InternalRational::divsame( InternalCF * c )
{
    return dividesame( c );
}

InternalCF* InternalRational::modulosame( InternalCF * c )
{
    return modsame( c );
}

InternalCF* InternalRational::modsame( InternalCF * )
{
    if ( deleteObject() ) delete this;
    return CFFactory::basic( 0 );
}

void InternalRational::divremsame( InternalCF * c, InternalCF*& quot, InternalCF*& rem )
{
    quot = copyObject();
    quot = quot->dividesame( c );
    rem = CFFactory::basic( 0 );
}

bool InternalRational::divremsamet( InternalCF* c, InternalCF*& quot, InternalCF*& rem )
{
    divremsame( c, quot, rem );
    return true;
}

//{{{ int InternalRational::comparesame, comparecoeff ( InternalCF * c )
//{{{ docu
//
// comparesame(), comparecoeff() - compare with an
//   InternalRational.
//
// comparesame() compares the CO=a/b and c=p/q using the
// equivalence a/b < p/q iff a*q < p*b.
//
// comparecoeff() compares the CO=a/b and the integer c using the
// equivalence a/b < c iff a < c*b.
//
// Note: Both operations may be relatively expensive due to the
// multiplications.
//
// See also: CanonicalForm::operator <(), CanonicalForm::operator ==()
//
//}}}
int
InternalRational::comparesame ( InternalCF * c )
{
    ASSERT( ! ::is_imm( c ) && c->levelcoeff() == RationalDomain, "incompatible base coefficients" );
    mpz_t dummy1, dummy2;
    mpz_init( dummy1 ); mpz_init( dummy2 );
    mpz_mul( dummy1, _num, MPQDEN( c ) );
    mpz_mul( dummy2, _den, MPQNUM( c ) );
    int result = mpz_cmp( dummy1, dummy2 );
    mpz_clear( dummy1 ); mpz_clear( dummy2 );
    return result;
}

int
InternalRational::comparecoeff ( InternalCF* c )
{
    if ( ::is_imm( c ) )
    {
        ASSERT( ::is_imm( c ) == INTMARK, "incompatible base coefficients" );
        mpz_t dummy;
        mpz_init_set_si( dummy, imm2int( c ) );
        mpz_mul( dummy, dummy, _den );
        int result = mpz_cmp( _num, dummy );
        mpz_clear( dummy );
        return result;
    }
    else
    {
        ASSERT( c->levelcoeff() == IntegerDomain, "incompatible base coefficients" );
        mpz_t dummy;
        mpz_init( dummy );
        mpz_mul( dummy, _den, InternalInteger::MPI( c ) );
        int result = mpz_cmp( _num, dummy );
        mpz_clear( dummy );
        return result;
    }
}
//}}}

InternalCF* InternalRational::addcoeff( InternalCF* c )
{
    ASSERT( ::is_imm( c ) == INTMARK || ! ::is_imm( c ), "expected integer" );
    mpz_t n, d;
    if ( ::is_imm( c ) )
    {
        int cc = imm2int( c );
        if ( cc == 0 )
            return this;
        else
        {
          mpz_init( n );
          if ( cc < 0 )
          {
            mpz_mul_ui( n, _den, -cc );
            mpz_sub( n, _num, n );
          }
          else
          {
            mpz_mul_ui( n, _den, cc );
            mpz_add( n, _num, n );
          }
        }
    }
    else
    {
        ASSERT( c->levelcoeff() == IntegerDomain, "expected integer" );
        mpz_init( n );
        mpz_mul( n, _den, InternalInteger::MPI( c ) );
        mpz_add( n, _num, n );
    }
    mpz_init_set( d, _den );
    // at this point there is no way that the result is not a true rational
    if ( deleteObject() ) delete this;
    return new InternalRational( n, d );
}

InternalCF* InternalRational::subcoeff( InternalCF* c, bool negate )
{
    ASSERT( ::is_imm( c ) == INTMARK || ! ::is_imm( c ), "expected integer" );
    mpz_t n, d;
    if ( ::is_imm( c ) )
    {
        int cc = imm2int( c );
        if ( cc == 0 )
        {
            if ( negate )
            {
                if ( getRefCount() == 1 )
                {
                    mpz_neg( _num, _num );
                    return this;
                }
                else
                {
                    decRefCount();
                    mpz_init_set( d, _den );
                    mpz_init_set( n, _num );
                    mpz_neg( n, n );
                    return new InternalRational( n, d );
                }
            }
            else
                return this;
        }
        mpz_init( n );
        if ( cc < 0 )
        {
            mpz_mul_ui( n, _den, -cc );
            mpz_neg( n, n );
        }
        else
            mpz_mul_ui( n, _den, cc );
        if ( negate )
            mpz_sub( n, n, _num );
        else
            mpz_sub( n, _num, n );
    }
    else
    {
        ASSERT( c->levelcoeff() == IntegerDomain, "expected integer" );
        mpz_init( n );
        mpz_mul( n, _den, InternalInteger::MPI( c ) );
        if ( negate )
            mpz_sub( n, n, _num );
        else
            mpz_sub( n, _num, n );
    }
    mpz_init_set( d, _den );
    // at this point there is no way that the result is not a true rational
    if ( deleteObject() ) delete this;
    return new InternalRational( n, d );
}

InternalCF* InternalRational::mulcoeff( InternalCF* c )
{
    ASSERT( ::is_imm( c ) == INTMARK || ! ::is_imm( c ), "expected integer" );
    mpz_t n, d, g;
    if ( ::is_imm( c ) )
    {
        int cc = imm2int( c );
        if ( cc == 0 )
        {
            if ( deleteObject() ) delete this;
            return CFFactory::basic( 0 );
        }
        mpz_init_set_si( n, cc );
    }
    else
    {
        ASSERT( c->levelcoeff() == IntegerDomain, "expected integer" );
        mpz_init_set( n, InternalInteger::MPI( c ) );
    }
    mpz_init( g );
    mpz_gcd( g, n, _den );
    if ( mpz_cmp_si( g, 1 ) == 0 )
    {
        mpz_mul( n, n, _num );
        mpz_init_set( d, _den );
    }
    else
    {
        mpz_divexact( n, n, g );
        mpz_mul( n, n, _num );
        mpz_init( d );
        mpz_divexact( d, _den, g );
    }
    mpz_clear( g );
    if ( deleteObject() ) delete this;
    if ( mpz_cmp_si( d, 1 ) == 0 )
    {
        mpz_clear( d );
        if ( mpz_is_imm( n ) )
        {
            InternalCF * res = int2imm( mpz_get_si( n ) );
            mpz_clear( n );
            return res;
        }
        else
        {
            return new InternalInteger( n );
        }
    }
    else
        return new InternalRational( n, d );
}

InternalCF* InternalRational::dividecoeff( InternalCF* c, bool invert )
{
    ASSERT( ::is_imm( c ) == INTMARK || ! ::is_imm( c ), "expected integer" );
    mpz_t n, d, g;
    if ( ::is_imm( c ) )
    {
        int cc = imm2int( c );
        ASSERT( c != 0 || invert, "divide by zero" );
        if ( cc == 0 )
        {
            // => invert
            if ( deleteObject() ) delete this;
            return CFFactory::basic( 0 );
        }
        if ( invert )
        {
            mpz_init_set_si( n, cc );
            mpz_mul( n, n, _den );
            mpz_init_set( d, _num );
        }
        else
        {
            mpz_init_set_si( d, cc );
            mpz_mul( d, d, _den );
            mpz_init_set( n, _num );
        }
    }
    else
    {
        ASSERT( c->levelcoeff() == IntegerDomain, "expected integer" );
        if ( invert )
        {
            mpz_init_set( n, InternalInteger::MPI( c ) );
            mpz_mul( n, n, _den );
            mpz_init_set( d, _num );
        }
        else
        {
            mpz_init_set( d, InternalInteger::MPI( c ) );
            mpz_mul( d, d, _den );
            mpz_init_set( n, _num );
        }
    }
    if ( mpz_cmp_si( d, 0 ) < 0 )
    {
        mpz_neg( d, d );
        mpz_neg( n, n );
    }
    mpz_init( g );
    mpz_gcd( g, n, d );
    if ( mpz_cmp_si( g, 1 ) != 0 )
    {
        mpz_divexact( d, d, g );
        mpz_divexact( n, n, g );
    }
    mpz_clear( g );
    if ( deleteObject() ) delete this;
    if ( ! invert )
    {
        // then there was no way for the result to become an integer
        return new InternalRational( n, d );
    }
    if ( mpz_cmp_si( d, 1 ) == 0 )
    {
        mpz_clear( d );
        if ( mpz_is_imm( n ) )
        {
            InternalCF * res = int2imm( mpz_get_si( n ) );
            mpz_clear( n );
            return res;
        }
        else
        {
            return new InternalInteger( n );
        }
    }
    else
        return new InternalRational( n, d );
}

InternalCF* InternalRational::divcoeff( InternalCF* c, bool invert )
{
    return dividecoeff( c, invert );
}

InternalCF* InternalRational::modulocoeff( InternalCF* c, bool invert )
{
    return modcoeff( c, invert );
}

InternalCF* InternalRational::modcoeff( InternalCF* c, bool invert )
{
    ASSERT( ::is_imm( c ) == INTMARK || ! ::is_imm( c ), "integer expected" );
    ASSERT( invert || ! ::is_imm( c ) || imm2int( c ) != 0, "divide by zero" );
    if ( deleteObject() ) delete this;
    return CFFactory::basic( 0 );
}

void InternalRational::divremcoeff( InternalCF* c, InternalCF*& quot, InternalCF*& rem, bool invert )
{
    quot = copyObject();
    quot = quot->dividecoeff( c, invert );
    rem = CFFactory::basic( 0 );
}

bool InternalRational::divremcoefft( InternalCF* c, InternalCF*& quot, InternalCF*& rem, bool invert )
{
    divremcoeff( c, quot, rem, invert );
    return true;
}

//{{{ InternalCF * InternalRational::bgcdsame, bgcdcoeff ( const InternalCF * const )
// docu: see CanonicalForm::bgcd()
InternalCF *
InternalRational::bgcdsame ( const InternalCF * const ) const
{
    return int2imm( 1 );
}

InternalCF *
InternalRational::bgcdcoeff ( const InternalCF * const )
{
    return int2imm( 1 );
}
//}}}

//{{{ InternalCF * InternalRational::bextgcdsame ( InternalCF * c, CanonicalForm & a, CanonicalForm & b )
// docu: see CanonicalForm::bextgcd()
InternalCF *
InternalRational::bextgcdsame ( InternalCF *, CanonicalForm & a, CanonicalForm & b )
{
    a = 1/CanonicalForm( copyObject() ); b = 0;
    return int2imm( 1 );
}

InternalCF *
InternalRational::bextgcdcoeff ( InternalCF *, CanonicalForm & a, CanonicalForm & b )
{
    a = 1/CanonicalForm( copyObject() ); b = 0;
    return int2imm( 1 );
}
//}}}

InternalCF * InternalRational::normalize_myself()
{
    ASSERT( getRefCount() == 1, "illegal operation" );
    mpz_t g;
    mpz_init( g );
    mpz_gcd( g, _num, _den );
    if ( mpz_cmp_si( g, 1 ) != 0 )
    {
        mpz_divexact( _num, _num, g );
        mpz_divexact( _den, _den, g );
    }
    // Hier brauchen wir ein mpz_clear, J.M.
    mpz_clear( g );
    if ( mpz_cmp_si( _den, 0 ) < 0 )
    {
        mpz_neg( _num, _num );
        mpz_neg( _den, _den );
    }
    if ( mpz_cmp_si( _den, 1 ) == 0 )
    {
        if ( mpz_is_imm( _num ) )
        {
            InternalCF * res = int2imm( mpz_get_si( _num ) );
            delete this;
            return res;
        }
        else
        {
            mpz_t res;
            mpz_init_set( res, _num );
            delete this;
            return new InternalInteger( res );
        }
    }
    else
        return this;
}


int InternalRational::intval() const
{
    ASSERT( mpz_cmp_si( _den, 1 ) == 0, "illegal operation" );
    return (int)mpz_get_si( _num );
}

//{{{ int InternalRational::sign () const
// docu: see CanonicalForm::sign()
int
InternalRational::sign () const
{
    return mpz_sgn( _num );
}
//}}}