-rw-r--r-- 5343 cryptattacktester-20230614/bigfloat_mpfi.cpp raw
#include <cassert>
#include <map>
#include <mpfr.h>
#include "bigfloat_mpfi.h"

using namespace std;

void bigfloat::set_default_precision(bigint p)
{
  mpfr_set_default_prec(p);
}

bigint bigfloat::get_default_precision(void)
{
  return mpfr_get_default_prec();
}

ostream &operator<<(ostream &o,const bigfloat &a)
{
  char *s = 0;
  assert(mpfr_asprintf(&s,"[%.9RDg,%.9RUg]",&a.x->left,&a.x->right) >= 0);
  o << s;
  mpfr_free_str(s);
  return o;
}

bigfloat bigfloat_guarantee_nonnegative(const bigfloat &a)
{
  bigfloat result;
  mpfr_t left,right;
  mpfr_init(left);
  mpfr_init(right);
  mpfi_get_left(left,a.x);
  mpfi_get_right(right,a.x);
  if (mpfr_sgn(left) < 0)
    mpfr_set_zero(left,1);
  // XXX: can produce empty interval
  mpfi_interv_fr(result.x,left,right);
  mpfr_clear(right);
  mpfr_clear(left);
  return result;
}

bool bigfloat::definitely_zero() { return mpfi_is_zero(x); }
bool bigfloat::definitely_positive() { return mpfi_is_strictly_pos(x); }

bool bigfloat::definitely_less(const bigfloat &b)
{
  // mpfi documentation says mpfi_cmp returns 1 for invalid operands
  return mpfi_cmp(x,b.x) < 0;
}

bool bigfloat::maybe_between(double f,double g)
{
  return mpfi_cmp_d(x,f) >= 0 && mpfi_cmp_d(x,g) <= 0;
}

// returns 1 if all of the following occur:
// * interval is a bounded interval [left,right]
// * left > 0
// * 2^reasonable(right-left) < left; i.e. right/left < 1+1/2^reasonable
bool bigfloat::definitely_positive_and_reasonably_precise(long reasonable)
{
  if (!mpfi_bounded_p(x)) return 0;
  if (!mpfi_is_strictly_pos(x)) return 0;
  mpfr_t left,gap;
  mpfr_init(left);
  mpfr_init(gap);
  mpfi_get_left(left,x);
  mpfi_get_right(gap,x);
  mpfr_sub(gap,gap,left,GMP_RNDU);
  mpfr_mul_2ui(gap,gap,reasonable,GMP_RNDU);
  bool result = mpfr_cmp(gap,left) < 0;
  mpfr_clear(gap);
  mpfr_clear(left);
  return result;
}

// returns total ordering on (left,right), suitable for sorting
// beware that NaN follows its own rules
bool operator<(const bigfloat &a,const bigfloat &b)
{
  mpfr_t aleft,aright,bleft,bright;
  mpfr_init(aleft);
  mpfr_init(aright);
  mpfr_init(bleft);
  mpfr_init(bright);
  mpfi_get_left(aleft,a.x);
  mpfi_get_right(aright,a.x);
  mpfi_get_left(bleft,b.x);
  mpfi_get_right(bright,b.x);
  bool result = mpfr_less_p(aleft,bleft) || (mpfr_equal_p(aleft,bleft) && mpfr_less_p(aright,bright));
  mpfr_clear(bright);
  mpfr_clear(bleft);
  mpfr_clear(aright);
  mpfr_clear(aleft);
  return result;
}

// returns total ordering on (left,right), suitable for sorting
// beware that NaN follows its own rules
bool operator>(const bigfloat &a,const bigfloat &b)
{
  mpfr_t aleft,aright,bleft,bright;
  mpfr_init(aleft);
  mpfr_init(aright);
  mpfr_init(bleft);
  mpfr_init(bright);
  mpfi_get_left(aleft,a.x);
  mpfi_get_right(aright,a.x);
  mpfi_get_left(bleft,b.x);
  mpfi_get_right(bright,b.x);
  bool result = mpfr_greater_p(aleft,bleft) || (mpfr_equal_p(aleft,bleft) && mpfr_greater_p(aright,bright));
  mpfr_clear(bright);
  mpfr_clear(bleft);
  mpfr_clear(aright);
  mpfr_clear(aleft);
  return result;
}

bigfloat operator+(const bigfloat &a,const bigfloat &b)
{
  bigfloat result;
  mpfi_add(result.x,a.x,b.x);
  return result;
}

bigfloat& operator+=(bigfloat &a,const bigfloat &b)
{
  mpfi_add(a.x,a.x,b.x);
  return a;
}

bigfloat operator-(const bigfloat &a,const bigfloat &b)
{
  bigfloat result;
  mpfi_sub(result.x,a.x,b.x);
  return result;
}

bigfloat operator-(const bigfloat &a)
{
  bigfloat result;
  mpfi_neg(result.x,a.x);
  return result;
}

bigfloat operator*(const bigfloat &a,const bigfloat &b)
{
  bigfloat result;
  mpfi_mul(result.x,a.x,b.x);
  return result;
}

bigfloat& operator*=(bigfloat &a,const bigfloat &b)
{
  mpfi_mul(a.x,a.x,b.x);
  return a;
}

bigfloat operator/(const bigfloat &a,const bigfloat &b)
{
  bigfloat result;
  mpfi_div(result.x,a.x,b.x);
  return result;
}

bigfloat log(const bigfloat &a) { bigfloat result; mpfi_log(result.x,a.x); return result; }
bigfloat log2(const bigfloat &a) { bigfloat result; mpfi_log2(result.x,a.x); return result; }
bigfloat log1p(const bigfloat &a) { bigfloat result; mpfi_log1p(result.x,a.x); return result; }
bigfloat exp(const bigfloat &a) { bigfloat result; mpfi_exp(result.x,a.x); return result; }
bigfloat expm1(const bigfloat &a) { bigfloat result; mpfi_expm1(result.x,a.x); return result; }
bigfloat exp2(const bigfloat &a) { bigfloat result; mpfi_exp2(result.x,a.x); return result; }

bigfloat pow(const bigfloat &a,const bigfloat &b)
{
  return exp(b*log(a));
}

bigint floor_as_bigint(const bigfloat &a)
{
  assert(mpfi_bounded_p(a.x)); // XXX
  mpz_t zleft;
  mpfr_t left;
  mpz_init(zleft);
  mpfr_init(left);
  mpfi_get_left(left,a.x);
  mpfr_get_z(zleft,left,GMP_RNDD);
  bigint result(zleft);
  mpfr_clear(left);
  mpz_clear(zleft);
  return result;
}

bigint ceil_as_bigint(const bigfloat &a)
{
  assert(mpfi_bounded_p(a.x)); // XXX
  mpz_t zright;
  mpfr_t right;
  mpz_init(zright);
  mpfr_init(right);
  mpfi_get_right(right,a.x);
  mpfr_get_z(zright,right,GMP_RNDU);
  bigint result(zright);
  mpfr_clear(right);
  mpz_clear(zright);
  return result;
}

bigfloat binomial_bigfloat(const bigfloat &n,const bigint &k)
{
  if (k < 0) return 0;
  if (k == 0) return 1;
  bigfloat n1 = n-1;
  bigint k1 = k-1;
  bigfloat result = binomial_bigfloat(n1,k1)*n/bigfloat(k);
  return result;
}