From 293dadb003aa547958c46d7f70d5b856b01f9f31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADtor=20Santos=20Costa?= Date: Fri, 28 May 2010 09:53:56 +0100 Subject: [PATCH] support for rational numbers make floor and friends return an integer (make it closer to SICStus). --- C/arith1.c | 179 +++++++------------------ C/bignum.c | 84 +++++++++++- C/globals.c | 3 +- C/gmp_support.c | 346 +++++++++++++++++++++++++++++++++++++++++++++--- C/write.c | 75 ++--------- H/Yapproto.h | 5 +- H/eval.h | 16 +++ H/iatoms.h | 2 + H/ratoms.h | 2 + H/tatoms.h | 4 + misc/ATOMS | 2 + 11 files changed, 495 insertions(+), 223 deletions(-) diff --git a/C/arith1.c b/C/arith1.c index fa2d2cf21..dfedb3784 100755 --- a/C/arith1.c +++ b/C/arith1.c @@ -37,9 +37,7 @@ float_to_int(Float v) if (i-v == 0.0) { return MkIntegerTerm(i); } else { - MP_INT o; - mpz_init_set_d(&o, v); - return Yap_MkBigIntTerm(&o); + return Yap_gmp_float_to_big(v); } #else return MkIntegerTerm(v); @@ -48,24 +46,6 @@ float_to_int(Float v) #define RBIG_FL(v) return(float_to_int(v)) -#if USE_GMP -static Term -process_iso_error(MP_INT *big, Term t, char *operation) -{ /* iso */ - Int sz = 2+mpz_sizeinbase(big,10); - char *s = Yap_AllocCodeSpace(sz); - - if (s != NULL) { - mpz_get_str(s, 10, big); - Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is %s(%s)", operation, s); - Yap_FreeCodeSpace(s); - RERROR(); - } else { - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is %s(t)",operation); - } -} -#endif - typedef struct init_un_eval { char *OpName; arith1_op f; @@ -104,7 +84,7 @@ get_float(Term t) { } #ifdef USE_GMP if (IsBigIntTerm(t)) { - return mpz_get_d(Yap_BigIntOfTerm(t)); + return Yap_gmp_to_float(t); } #endif return 0.0; @@ -434,33 +414,13 @@ eval1(Int fi, Term t) { switch (ETypeOfTerm(t)) { case long_int_e: - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is floor(%f)", IntegerOfTerm(t)); - } else { - RFLOAT(IntegerOfTerm(t)); - } + return t; case double_e: dbl = FloatOfTerm(t); break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - MP_INT *big = Yap_BigIntOfTerm(t); - Int sz = 2+mpz_sizeinbase(big,10); - char *s = Yap_AllocCodeSpace(sz); - - if (s != NULL) { - mpz_get_str(s, 10, big); - Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is floor(%s)", s); - Yap_FreeCodeSpace(s); - RERROR(); - } else { - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is floor(t)"); - } - } else { - dbl = mpz_get_d(Yap_BigIntOfTerm(t)); - } - break; + return Yap_gmp_floor(t); #endif default: RERROR(); @@ -483,22 +443,13 @@ eval1(Int fi, Term t) { Float dbl; switch (ETypeOfTerm(t)) { case long_int_e: - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is ceiling(%f)", IntegerOfTerm(t)); - } else { - RFLOAT(IntegerOfTerm(t)); - } + return t; case double_e: dbl = FloatOfTerm(t); break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return process_iso_error(Yap_BigIntOfTerm(t), t, "ceiling"); - } else { - dbl = mpz_get_d(Yap_BigIntOfTerm(t)); - } - break; + return Yap_gmp_ceiling(t); #endif default: RERROR(); @@ -522,21 +473,13 @@ eval1(Int fi, Term t) { switch (ETypeOfTerm(t)) { case long_int_e: - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is round(%ld)", IntegerOfTerm(t)); - } else { - return t; - } + return t; case double_e: dbl = FloatOfTerm(t); break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { - return process_iso_error(Yap_BigIntOfTerm(t), t, "round"); - } - return t; - break; + return Yap_gmp_round(t); #endif default: RERROR(); @@ -560,19 +503,13 @@ eval1(Int fi, Term t) { Float dbl; switch (ETypeOfTerm(t)) { case long_int_e: - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is round(%ld)", IntegerOfTerm(t)); - } return t; case double_e: dbl = FloatOfTerm(t); break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is round(BIGNUM)"); - } - return t; + return Yap_gmp_trunc(t); #endif default: RERROR(); @@ -588,18 +525,10 @@ eval1(Int fi, Term t) { (%f)",dbl); } #endif - if (dbl <= (Float)Int_MAX && dbl >= (Float)Int_MIN) { - RINT((Int) dbl); - } else { -#ifdef USE_GMP - MP_INT new; - - mpz_init_set_d(&new, dbl); - RBIG(&new); -#else - return Yap_ArithError(EVALUATION_ERROR_INT_OVERFLOW, MkFloatTerm(dbl), "integer/1"); -#endif - } + if (dbl < 0.0) + RBIG_FL(ceil(dbl)); + else + RBIG_FL(floor(dbl)); } case op_float: switch (ETypeOfTerm(t)) { @@ -609,11 +538,37 @@ eval1(Int fi, Term t) { return t; case big_int_e: #ifdef USE_GMP - RFLOAT(mpz_get_d(Yap_BigIntOfTerm(t))); + RFLOAT(Yap_gmp_to_float(t)); #endif default: RERROR(); } + case op_rational: + switch (ETypeOfTerm(t)) { + case long_int_e: + return t; +#ifdef USE_GMP + case double_e: + return Yap_gmp_float_to_rational(FloatOfTerm(t)); +#endif + case big_int_e: + return t; + default: + RERROR(); + } + case op_rationalize: + switch (ETypeOfTerm(t)) { + case long_int_e: + return t; +#ifdef USE_GMP + case double_e: + return Yap_gmp_float_rationalize(FloatOfTerm(t)); +#endif + case big_int_e: + return t; + default: + RERROR(); + } case op_abs: switch (ETypeOfTerm(t)) { case long_int_e: @@ -622,13 +577,7 @@ eval1(Int fi, Term t) { RFLOAT(fabs(FloatOfTerm(t))); case big_int_e: #ifdef USE_GMP - { - MP_INT new; - - mpz_init_set(&new, Yap_BigIntOfTerm(t)); - mpz_abs(&new, &new); - RBIG(&new); - } + return Yap_gmp_abs_big(t); #endif default: RERROR(); @@ -641,13 +590,7 @@ eval1(Int fi, Term t) { return Yap_ArithError(TYPE_ERROR_INTEGER, t, "msb(%f)", FloatOfTerm(t)); case big_int_e: #ifdef USE_GMP - { MP_INT *big = Yap_BigIntOfTerm(t); - if ( mpz_sgn(big) <= 0 ) { - return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, - "msb/1 received bignum"); - } - RINT(mpz_sizeinbase(big,2)); - } + return Yap_gmp_msb(t); #endif default: RERROR(); @@ -660,13 +603,7 @@ eval1(Int fi, Term t) { return Yap_ArithError(TYPE_ERROR_INTEGER, t, "lsb(%f)", FloatOfTerm(t)); case big_int_e: #ifdef USE_GMP - { MP_INT *big = Yap_BigIntOfTerm(t); - if ( mpz_sgn(big) <= 0 ) { - return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, - "lsb/1 received bignum"); - } - RINT(mpz_scan1(big,0)); - } + return Yap_gmp_lsb(t); #endif default: RERROR(); @@ -679,13 +616,7 @@ eval1(Int fi, Term t) { return Yap_ArithError(TYPE_ERROR_INTEGER, t, "popcount(%f)", FloatOfTerm(t)); case big_int_e: #ifdef USE_GMP - { MP_INT *big = Yap_BigIntOfTerm(t); - if ( mpz_sgn(big) <= 0 ) { - return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, - "popcount/1 received negative bignum"); - } - RINT(mpz_popcount(big)); - } + return Yap_gmp_popcount(t); #endif default: RERROR(); @@ -707,11 +638,7 @@ eval1(Int fi, Term t) { break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return process_iso_error(Yap_BigIntOfTerm(t), t, "float_fractional_part"); - } else { - RFLOAT(0.0); - } + return Yap_gmp_float_fractional_part(t); #endif default: RERROR(); @@ -719,21 +646,13 @@ eval1(Int fi, Term t) { case op_fintp: switch (ETypeOfTerm(t)) { case long_int_e: - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is float_integer_part(%f)", IntegerOfTerm(t)); - } else { - RFLOAT(IntegerOfTerm(t)); - } + return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is float_integer_part(%f)", IntegerOfTerm(t)); case double_e: RFLOAT(rint(FloatOfTerm(t))); break; case big_int_e: #ifdef USE_GMP - if (yap_flags[LANGUAGE_MODE_FLAG] == 1) { /* iso */ - return process_iso_error(Yap_BigIntOfTerm(t), t, "float_integer_part"); - } else { - RFLOAT(mpz_get_d(Yap_BigIntOfTerm(t))); - } + return Yap_gmp_float_integer_part(t); #endif default: RERROR(); @@ -755,7 +674,7 @@ eval1(Int fi, Term t) { } case big_int_e: #ifdef USE_GMP - RINT(mpz_sgn(Yap_BigIntOfTerm(t))); + return Yap_gmp_sign(t); #endif default: RERROR(); @@ -818,6 +737,8 @@ static InitUnEntry InitUnTab[] = { {"lgamma", op_lgamma}, {"erf",op_erf}, {"erfc",op_erfc}, + {"rational",op_rational}, + {"rationalize",op_rationalize}, {"random", op_random1} }; diff --git a/C/bignum.c b/C/bignum.c index 1b65e23db..3430c74bc 100755 --- a/C/bignum.c +++ b/C/bignum.c @@ -83,7 +83,7 @@ Yap_MkBigRatTerm(MP_RAT *big) } H[0] = (CELL)FunctorBigInt; H[1] = BIG_RATIONAL; - dst->_mp_alloc = 0; + dst->_mp_size = 0; rat = (MP_RAT *)(dst+1); rat->_mp_num._mp_size = num->_mp_size; rat->_mp_num._mp_alloc = num->_mp_alloc; @@ -95,7 +95,7 @@ Yap_MkBigRatTerm(MP_RAT *big) nlimbs = (den->_mp_alloc)*(sizeof(mp_limb_t)/CellSize); memmove((void *)(H), (const void *)(den->_mp_d), nlimbs*CellSize); H += nlimbs; - dst->_mp_size = (H-(CELL *)rat); + dst->_mp_alloc = (H-(CELL *)(dst+1)); H[0] = EndSpecials; H++; return AbsAppl(ret); @@ -113,6 +113,17 @@ Yap_BigRatOfTerm(Term t) return new; } +Term +Yap_RatTermToApplTerm(Term t) +{ + Term ts[2]; + MP_RAT *rat = Yap_BigRatOfTerm(t); + + ts[0] = Yap_MkBigIntTerm(mpq_numref(rat)); + ts[1] = Yap_MkBigIntTerm(mpq_denref(rat)); + return Yap_MkApplTerm(FunctorRDiv,2,ts); +} + #endif @@ -150,7 +161,11 @@ p_is_bignum(void) #ifdef USE_GMP Term t = Deref(ARG1); return( - IsNonVarTerm(t) && IsApplTerm(t) && FunctorOfTerm(t) == FunctorBigInt); + IsNonVarTerm(t) && + IsApplTerm(t) && + FunctorOfTerm(t) == FunctorBigInt && + RepAppl(t)[1] == BIG_INT + ); #else return FALSE; #endif @@ -166,9 +181,72 @@ p_has_bignums(void) #endif } +static Int +p_is_rational(void) +{ + Term t = Deref(ARG1); + if (IsVarTerm(t)) + return FALSE; + if (IsIntTerm(t)) + return TRUE; + if (IsApplTerm(t)) { + Functor f = FunctorOfTerm(t); + CELL *pt; + + if (f == FunctorLongInt) + return TRUE; + if (f != FunctorBigInt) + return FALSE; + pt = RepAppl(t); + return ( pt[1] == BIG_RATIONAL || pt[1] == BIG_INT ); + } + return FALSE; +} + +static Int +p_rational(void) +{ +#ifdef USE_GMP + Term t = Deref(ARG1); + Functor f; + CELL *pt; + MP_RAT *rat; + Term t1, t2; + + if (IsVarTerm(t)) + return FALSE; + if (!IsApplTerm(t)) + return FALSE; + f = FunctorOfTerm(t); + if (f != FunctorBigInt) + return FALSE; + pt = RepAppl(t); + if (pt[1] != BIG_RATIONAL) + return FALSE; + rat = Yap_BigRatOfTerm(t); + while ((t1 = Yap_MkBigIntTerm(mpq_numref(rat))) == TermNil || + (t2 = Yap_MkBigIntTerm(mpq_denref(rat))) == TermNil) { + UInt size = + (mpq_numref(rat)->_mp_alloc)*(sizeof(mp_limb_t)/CellSize) + + (mpq_denref(rat)->_mp_alloc)*(sizeof(mp_limb_t)/CellSize); + if (!Yap_gcl(size, 3, ENV, P)) { + Yap_Error(OUT_OF_STACK_ERROR, t, Yap_ErrorMessage); + return FALSE; + } + } + return + Yap_unify(ARG2, t1) && + Yap_unify(ARG3, t2); +#else + return FALSE; +#endif +} + void Yap_InitBigNums(void) { Yap_InitCPred("$has_bignums", 0, p_has_bignums, SafePredFlag|HiddenPredFlag); Yap_InitCPred("$bignum", 1, p_is_bignum, SafePredFlag|HiddenPredFlag); + Yap_InitCPred("rational", 3, p_rational, 0); + Yap_InitCPred("rational", 1, p_is_rational, SafePredFlag); } diff --git a/C/globals.c b/C/globals.c index 4561baf05..dcc9626d3 100644 --- a/C/globals.c +++ b/C/globals.c @@ -398,7 +398,8 @@ copy_complex_term(register CELL *pt0, register CELL *pt0_end, int share, int cop default: { /* big int */ - UInt sz = ArenaSz(d0), i; + UInt sz = (sizeof(MP_INT)+3*CellSize+ + ((MP_INT *)(ap2+2))->_mp_alloc*sizeof(mp_limb_t))/CellSize, i; if (H > ASP - (MIN_ARENA_SIZE+sz)) { goto overflow; diff --git a/C/gmp_support.c b/C/gmp_support.c index c0704e0e0..51e894ea2 100755 --- a/C/gmp_support.c +++ b/C/gmp_support.c @@ -462,10 +462,13 @@ Yap_gmp_mul_big_big(Term t1, Term t2) MP_RAT new; MP_RAT *b1, bb1; MP_RAT *b2, bb2; + int f1 = FALSE, f2 = FALSE; + if (pt1[1] == BIG_INT) { b1 = &bb1; mpq_init(b1); mpq_set_z(b1, Yap_BigIntOfTerm(t1)); + f1 = TRUE; } else { b1 = Yap_BigRatOfTerm(t1); } @@ -473,11 +476,14 @@ Yap_gmp_mul_big_big(Term t1, Term t2) b2 = &bb2; mpq_init(b2); mpq_set_z(b2, Yap_BigIntOfTerm(t2)); + f2 = TRUE; } else { b2 = Yap_BigRatOfTerm(t2); } mpq_init(&new); mpq_mul(&new, b1, b2); + if (f1) mpq_clear(b1); + if (f2) mpq_clear(b2); return MkRatAndClose(&new); } } @@ -751,6 +757,15 @@ Yap_gmp_gcd_int_big(Int i, Term t) } } +Term +Yap_gmp_float_to_big(Float v) +{ + MP_INT new; + + mpz_init_set_d(&new, v); + return MkBigAndClose(&new); +} + Float Yap_gmp_to_float(Term t) { @@ -939,6 +954,10 @@ Yap_gmq_rdiv_int_int(Int i1, Int i2) MP_RAT new; mpq_init(&new); + if (i2 < 0) { + i1 = -i1; + i2 = -i2; + } mpq_set_si(&new, i1, i2); mpq_canonicalize(&new); return MkRatAndClose(&new); @@ -947,46 +966,77 @@ Yap_gmq_rdiv_int_int(Int i1, Int i2) Term Yap_gmq_rdiv_int_big(Int i1, Term t2) { - MP_RAT new, new2; - MP_INT *b = Yap_BigIntOfTerm(t2); - + MP_RAT new; + CELL *pt2 = RepAppl(t2); mpq_init(&new); mpq_set_si(&new, i1, 1L); - mpq_init(&new2); - mpq_set_z(&new2, b); - mpq_div(&new,&new,&new2); - mpq_clear(&new2); + if (pt2[1] == BIG_INT) { + MP_RAT new2; + MP_INT *b = Yap_BigIntOfTerm(t2); + + mpq_init(&new2); + mpq_set_z(&new2, b); + mpq_div(&new,&new,&new2); + mpq_clear(&new2); + } else { + MP_RAT *b = Yap_BigRatOfTerm(t2); + mpq_div(&new,&new,b); + } return MkRatAndClose(&new); } Term Yap_gmq_rdiv_big_int(Term t1, Int i2) { - MP_RAT new, new2; - MP_INT *b = Yap_BigIntOfTerm(t1); + MP_RAT new; + CELL *pt1 = RepAppl(t1); mpq_init(&new); mpq_set_si(&new, i2, 1L); - mpq_init(&new2); - mpq_set_z(&new2, b); - mpq_div(&new,&new2,&new); - mpq_clear(&new2); + if (pt1[1] == BIG_INT) { + MP_INT *b = Yap_BigIntOfTerm(t1); + MP_RAT new2; + + mpq_init(&new2); + mpq_set_z(&new2, b); + mpq_div(&new,&new2,&new); + mpq_clear(&new2); + } else { + MP_RAT *b = Yap_BigRatOfTerm(t1); + + mpq_div(&new,b,&new); + } return MkRatAndClose(&new); } Term Yap_gmq_rdiv_big_big(Term t1, Term t2) { - MP_RAT new, new2; - MP_INT *b1 = Yap_BigIntOfTerm(t1); - MP_INT *b2 = Yap_BigIntOfTerm(t2); + MP_RAT new; + CELL *pt1 = RepAppl(t1); + CELL *pt2 = RepAppl(t2); mpq_init(&new); - mpq_set_z(&new, b1); - mpq_init(&new2); - mpq_set_z(&new2, b2); - mpq_div(&new,&new,&new2); - mpq_clear(&new2); + if (pt1[1] == BIG_INT) { + MP_INT *b1 = Yap_BigIntOfTerm(t1); + mpq_set_z(&new, b1); + } else { + MP_RAT *b1 = Yap_BigRatOfTerm(t1); + mpq_set(&new, b1); + } + + if (pt2[1] == BIG_INT) { + MP_RAT new2; + MP_INT *b2 = Yap_BigIntOfTerm(t2); + + mpq_init(&new2); + mpq_set_z(&new2, b2); + mpq_div(&new,&new,&new2); + mpq_clear(&new2); + } else { + MP_RAT *b2 = Yap_BigRatOfTerm(t2); + mpq_div(&new,&new,b2); + } return MkRatAndClose(&new); } @@ -1157,13 +1207,82 @@ Yap_gmp_neg_big(Term t) return MkBigAndClose(&new); } else { MP_RAT *b = Yap_BigRatOfTerm(t); - MP_INT new; + MP_RAT new; mpq_init(&new); mpq_neg(&new, b); return MkRatAndClose(&new); } } +Term +Yap_gmp_float_to_rational(Float dbl) +{ + MP_RAT new; + mpq_init(&new); + mpq_set_d(&new, dbl); + return MkRatAndClose(&new); +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +A is rationalize(Float) + +Introduced on the suggestion of Richard O'Keefe after the Common Lisp +standard. The algorithm is taken from figure 3 in ``A Rational Rotation +Method for Robust Geometric Algorithms'' by John Canny, Bruce Donald and +Eugene K. Ressler. Found at + +http://www.cs.dartmouth.edu/~brd/papers/rotations-scg92.pdf +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +#ifndef DBL_EPSILON /* normal for IEEE 64-bit double */ +#define DBL_EPSILON 0.00000000000000022204 +#endif + +Term +Yap_gmp_float_rationalize(Float dbl) +{ + Float e0 = dbl, p0 = 0.0, q0 = 1.0; + Float e1 = -1.0, p1 = 1.0, q1 = 0.0; + Float d; + MP_RAT new; + + do { Float r = floor(e0/e1); + Float e00 = e0, p00 = p0, q00 = q0; + e0 = e1; + p0 = p1; + q0 = q1; + e1 = e00 - r*e1; + p1 = p00 - r*p1; + q1 = q00 - r*q1; + + d = p1/q1 - dbl; + } while(fabs(d) > DBL_EPSILON); + + mpz_init_set_d(mpq_numref(&new), p1); + mpz_init_set_d(mpq_denref(&new), q1); + mpq_canonicalize(&new); /* is this needed? */ + return MkRatAndClose(&new); +} + +Term +Yap_gmp_abs_big(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + MP_INT *b = Yap_BigIntOfTerm(t); + MP_INT new; + mpz_init_set(&new, b); + mpz_abs(&new, &new); + return MkBigAndClose(&new); + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_RAT new; + mpq_init(&new); + mpq_abs(&new, b); + return MkRatAndClose(&new); + } +} + Term Yap_gmp_unot_big(Term t) { @@ -1179,6 +1298,189 @@ Yap_gmp_unot_big(Term t) } } +Term +Yap_gmp_floor(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return t; + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_INT new; + mpz_init(&new); + mpz_set_q(&new, b); + if (mpq_sgn(b) < 0 && mpz_cmp_si(mpq_denref(b),1L) != 0) { + mpz_sub_ui(&new,&new,1L); + } + return MkBigAndClose(&new); + } +} + +Term +Yap_gmp_ceiling(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return t; + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_INT new; + mpz_init(&new); + mpz_set_q(&new, b); + if (mpq_sgn(b) > 0 && mpz_cmp_si(mpq_denref(b),1L) != 0) { + mpz_add_ui(&new,&new,1L); + } + return MkBigAndClose(&new); + } +} + +Term +Yap_gmp_round(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return t; + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_INT new; + MP_RAT half, q; + + mpq_init(&half); + mpq_init(&q); + mpq_set_ui(&half, 1, 2); /* 1/2 */ + if ( mpq_sgn(b) > 0 ) + mpq_add(&q, b, &half); + else { + mpq_sub(&q, b, &half); + } + mpz_init(&new); + mpz_set_q(&new, &q); + mpq_clear(&half); + mpq_clear(&q); + return MkBigAndClose(&new); + } +} + +Term +Yap_gmp_trunc(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return t; + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_INT new; + int sgn = mpq_sgn(b); + + if (sgn) + mpq_neg(b, b); + mpz_init(&new); + mpz_set_q(&new, b); + if (sgn) { + mpq_neg(b, b); + mpz_neg(&new, &new); + } + return MkBigAndClose(&new); + } +} + +Term +Yap_gmp_float_fractional_part(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is float_fractional_part(%f)", FloatOfTerm(t)); + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_RAT new; + + mpq_init(&new); + mpz_tdiv_q(mpq_numref(&new), + mpq_numref(b), + mpq_denref(b)); + mpz_set_ui(mpq_denref(&new), 1); + mpq_sub(&new, b, &new); + return MkRatAndClose(&new); + } +} + +Term +Yap_gmp_float_integer_part(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return Yap_ArithError(TYPE_ERROR_FLOAT, t, "X is float_integer_part(%f)", FloatOfTerm(t)); + } else { + MP_RAT *b = Yap_BigRatOfTerm(t); + MP_INT new; + + mpz_init(&new); + mpz_tdiv_q(&new, + mpq_numref(b), + mpq_denref(b)); + return MkBigAndClose(&new); + } +} + +Term +Yap_gmp_sign(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + return MkIntegerTerm(mpz_sgn(Yap_BigIntOfTerm(t))); + } else { + return MkIntegerTerm(mpq_sgn(Yap_BigRatOfTerm(t))); + } +} + +Term +Yap_gmp_lsb(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + MP_INT *big = Yap_BigIntOfTerm(t); + if ( mpz_sgn(big) <= 0 ) { + return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, + "lsb/1 received negative bignum"); + } + return MkIntegerTerm(mpz_scan1(big,0)); + } else { + return Yap_ArithError(TYPE_ERROR_INTEGER, t, "lsb"); + } +} + +Term +Yap_gmp_msb(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + MP_INT *big = Yap_BigIntOfTerm(t); + if ( mpz_sgn(big) <= 0 ) { + return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, + "msb/1 received negative bignum"); + } + return MkIntegerTerm(mpz_sizeinbase(big,2)); + } else { + return Yap_ArithError(TYPE_ERROR_INTEGER, t, "popcount"); + } +} + +Term +Yap_gmp_popcount(Term t) +{ + CELL *pt = RepAppl(t); + if (pt[1] == BIG_INT) { + MP_INT *big = Yap_BigIntOfTerm(t); + if ( mpz_sgn(big) <= 0 ) { + return Yap_ArithError(DOMAIN_ERROR_NOT_LESS_THAN_ZERO, t, + "popcount/1 received negative bignum"); + } + return MkIntegerTerm(mpz_popcount(big)); + } else { + return Yap_ArithError(TYPE_ERROR_INTEGER, t, "popcount"); + } +} + #endif diff --git a/C/write.c b/C/write.c index d2f0935d8..76602c7d9 100755 --- a/C/write.c +++ b/C/write.c @@ -191,7 +191,6 @@ write_mpq(MP_RAT *q, wrf writewch) { char *s; size_t sz; - fprintf(stderr,"%ld %ld\n",mpz_sizeinbase (mpq_numref(q), 10),mpz_sizeinbase (mpq_denref(q), 10)); sz = ((size_t)3) +mpz_sizeinbase(mpq_numref(q), 10)+ mpz_sizeinbase (mpq_denref(q), 10); s = ensure_space(sz); if (mpq_sgn(q) < 0) { @@ -214,23 +213,24 @@ write_mpq(MP_RAT *q, wrf writewch) { } #endif + /* writes a bignum */ static void -writebig(Term t, wrf writewch) /* writes an integer */ +writebig(Term t, int p, int depth, int rinfixarg, struct write_globs *wglb, struct rewind_term *rwt) { #ifdef USE_GMP CELL *pt = RepAppl(t)+1; if (pt[0] == BIG_INT) { MP_INT *big = Yap_BigIntOfTerm(t); - write_mpint(big, writewch); + write_mpint(big, wglb->writewch); return; } else if (pt[0] == BIG_RATIONAL) { - MP_RAT *q = Yap_BigRatOfTerm(t); - write_mpq(q, writewch); + Term trat = Yap_RatTermToApplTerm(t); + writeTerm(trat, p, depth, rinfixarg, wglb, rwt); return; } #endif - wrputs("0",writewch); + wrputs("0",wglb->writewch); } static void @@ -739,66 +739,9 @@ writeTerm(Term t, int p, int depth, int rinfixarg, struct write_globs *wglb, str case (CELL)FunctorLongInt: wrputn(LongIntOfTerm(t),wglb->writewch); return; - case (CELL)FunctorBigInt: - writebig(t,wglb->writewch); - return; -#ifdef USE_GMP - { - MP_INT *big = Yap_BigIntOfTerm(t); - char *s; - s = (char *) Yap_PreAllocCodeSpace(); - while (s+3+mpz_sizeinbase(big, 10) >= (char *)AuxSp) { -#if USE_SYSTEM_MALLOC - /* may require stack expansion */ - if (!Yap_ExpandPreAllocCodeSpace(3+mpz_sizeinbase(big, 10), NULL, TRUE)) { - s = NULL; - break; - } - s = (char *) Yap_PreAllocCodeSpace(); -#else - s = NULL; -#endif - } - if (!s) { - s = (char *)TR; - while (s+3+mpz_sizeinbase(big, 10) >= Yap_TrailTop) { - if (!Yap_growtrail((3+mpz_sizeinbase(big, 10))/sizeof(CELL), FALSE)) { - s = NULL; - break; - } - s = (char *)TR; - } - } - if (!s) { - s = (char *)H; - if (s+3+mpz_sizeinbase(big, 10) >= (char *)ASP) { - Yap_Error(OUT_OF_STACK_ERROR,TermNil,"not enough space to write bignum: it requires %d bytes", 3+mpz_sizeinbase(big, 10)); - s = NULL; - } - } - if (mpz_sgn(big) < 0) { - if (lastw == symbol) - wrputc(' ', wglb->writewch); - } else { - if (lastw == alphanum) - wrputc(' ', wglb->writewch); - } - if (!s) { - s = mpz_get_str(NULL, 10, big); - if (!s) - return; - wrputs(s,wglb->writewch); - free(s); - } else { - mpz_get_str(s, 10, big); - wrputs(s,wglb->writewch); - } - } -#else - { - wrputs("0",wglb->writewch); - } -#endif + /* case (CELL)FunctorBigInt: */ + default: + writebig(t, p, depth, rinfixarg, wglb, rwt); return; } } diff --git a/H/Yapproto.h b/H/Yapproto.h index 44c49cca2..576a17351 100755 --- a/H/Yapproto.h +++ b/H/Yapproto.h @@ -117,8 +117,9 @@ void STD_PROTO(Yap_InitAttVarPreds,(void)); void STD_PROTO(Yap_InitBBPreds,(void)); /* bignum.c */ -Term STD_PROTO(Yap_MkULLIntTerm,(YAP_ULONG_LONG)); -void STD_PROTO(Yap_InitBigNums,(void)); +Term STD_PROTO(Yap_MkULLIntTerm, (YAP_ULONG_LONG)); +Term STD_PROTO(Yap_RatTermToApplTerm, (Term)); +void STD_PROTO(Yap_InitBigNums, (void)); /* c_interface.c */ Int STD_PROTO(YAP_Execute,(struct pred_entry *, CPredicate)); diff --git a/H/eval.h b/H/eval.h index 3ce24a6ac..2f9564a8f 100644 --- a/H/eval.h +++ b/H/eval.h @@ -95,6 +95,8 @@ typedef enum { op_lgamma, op_erf, op_erfc, + op_rational, + op_rationalize, op_random1 } arith1_op; @@ -263,6 +265,9 @@ Term STD_PROTO(Yap_gmp_gcd_big_big,(Term,Term)); Term STD_PROTO(Yap_gmp_big_from_64bits,(YAP_LONG_LONG)); +Term STD_PROTO(Yap_gmp_float_to_big,(Float)); +Term STD_PROTO(Yap_gmp_float_to_rational,(Float)); +Term STD_PROTO(Yap_gmp_float_rationalize,(Float)); Float STD_PROTO(Yap_gmp_to_float,(Term)); Term STD_PROTO(Yap_gmp_add_float_big,(Float, Term)); Term STD_PROTO(Yap_gmp_sub_float_big,(Float, Term)); @@ -278,8 +283,19 @@ int STD_PROTO(Yap_gmp_cmp_big_float,(Term, Float)); int STD_PROTO(Yap_gmp_cmp_big_big,(Term, Term)); Term STD_PROTO(Yap_gmp_neg_int,(Int)); +Term STD_PROTO(Yap_gmp_abs_big,(Term)); Term STD_PROTO(Yap_gmp_neg_big,(Term)); Term STD_PROTO(Yap_gmp_unot_big,(Term)); +Term STD_PROTO(Yap_gmp_floor,(Term)); +Term STD_PROTO(Yap_gmp_ceiling,(Term)); +Term STD_PROTO(Yap_gmp_round,(Term)); +Term STD_PROTO(Yap_gmp_trunc,(Term)); +Term STD_PROTO(Yap_gmp_float_fractional_part,(Term)); +Term STD_PROTO(Yap_gmp_float_integer_part,(Term)); +Term STD_PROTO(Yap_gmp_sign,(Term)); +Term STD_PROTO(Yap_gmp_lsb,(Term)); +Term STD_PROTO(Yap_gmp_msb,(Term)); +Term STD_PROTO(Yap_gmp_popcount,(Term)); #endif diff --git a/H/iatoms.h b/H/iatoms.h index 3d186ea6b..e23025b03 100644 --- a/H/iatoms.h +++ b/H/iatoms.h @@ -226,6 +226,7 @@ AtomRepeatSpace = Yap_LookupAtom("repeat "); AtomReposition = Yap_LookupAtom("reposition"); AtomRepresentationError = Yap_LookupAtom("representation_error"); + AtomRDiv = Yap_LookupAtom("rdiv"); AtomResize = Yap_LookupAtom("resize"); AtomResourceError = Yap_LookupAtom("resource_error"); AtomRestoreRegs = Yap_FullLookupAtom("$restore_regs"); @@ -389,6 +390,7 @@ FunctorPrologConstraint = Yap_MkFunctor(AtomProlog,2); FunctorQuery = Yap_MkFunctor(AtomQuery,1); FunctorRecordedWithKey = Yap_MkFunctor(AtomRecordedWithKey,6); + FunctorRDiv = Yap_MkFunctor(AtomRDiv,2); FunctorRedoFreeze = Yap_MkFunctor(AtomRedoFreeze,3); FunctorRepresentationError = Yap_MkFunctor(AtomRepresentationError,1); FunctorResourceError = Yap_MkFunctor(AtomResourceError,1); diff --git a/H/ratoms.h b/H/ratoms.h index 3f4f568fc..a53f76c3c 100644 --- a/H/ratoms.h +++ b/H/ratoms.h @@ -226,6 +226,7 @@ AtomRepeatSpace = AtomAdjust(AtomRepeatSpace); AtomReposition = AtomAdjust(AtomReposition); AtomRepresentationError = AtomAdjust(AtomRepresentationError); + AtomRDiv = AtomAdjust(AtomRDiv); AtomResize = AtomAdjust(AtomResize); AtomResourceError = AtomAdjust(AtomResourceError); AtomRestoreRegs = AtomAdjust(AtomRestoreRegs); @@ -389,6 +390,7 @@ FunctorPrologConstraint = FuncAdjust(FunctorPrologConstraint); FunctorQuery = FuncAdjust(FunctorQuery); FunctorRecordedWithKey = FuncAdjust(FunctorRecordedWithKey); + FunctorRDiv = FuncAdjust(FunctorRDiv); FunctorRedoFreeze = FuncAdjust(FunctorRedoFreeze); FunctorRepresentationError = FuncAdjust(FunctorRepresentationError); FunctorResourceError = FuncAdjust(FunctorResourceError); diff --git a/H/tatoms.h b/H/tatoms.h index 00d3b3d0b..3710678af 100644 --- a/H/tatoms.h +++ b/H/tatoms.h @@ -450,6 +450,8 @@ #define AtomReposition Yap_heap_regs->AtomReposition_ Atom AtomRepresentationError_; #define AtomRepresentationError Yap_heap_regs->AtomRepresentationError_ + Atom AtomRDiv_; +#define AtomRDiv Yap_heap_regs->AtomRDiv_ Atom AtomResize_; #define AtomResize Yap_heap_regs->AtomResize_ Atom AtomResourceError_; @@ -776,6 +778,8 @@ #define FunctorQuery Yap_heap_regs->FunctorQuery_ Functor FunctorRecordedWithKey_; #define FunctorRecordedWithKey Yap_heap_regs->FunctorRecordedWithKey_ + Functor FunctorRDiv_; +#define FunctorRDiv Yap_heap_regs->FunctorRDiv_ Functor FunctorRedoFreeze_; #define FunctorRedoFreeze Yap_heap_regs->FunctorRedoFreeze_ Functor FunctorRepresentationError_; diff --git a/misc/ATOMS b/misc/ATOMS index d59578e7e..35d7e439a 100644 --- a/misc/ATOMS +++ b/misc/ATOMS @@ -231,6 +231,7 @@ A Repeat N "repeat" A RepeatSpace N "repeat " A Reposition N "reposition" A RepresentationError N "representation_error" +A RDiv N "rdiv" A Resize N "resize" A ResourceError N "resource_error" A RestoreRegs F "$restore_regs" @@ -394,6 +395,7 @@ F Portray Portray 1 F PrologConstraint Prolog 2 F Query Query 1 F RecordedWithKey RecordedWithKey 6 +F RDiv RDiv 2 F RedoFreeze RedoFreeze 3 F RepresentationError RepresentationError 1 F ResourceError ResourceError 1