libfplll-dev 5.2.0-3build1 / usr / include / fplll / nr / dpe.h

/* Copyright (C) 2004, 2005, 2006, 2008 Patrick Pelissier, Paul Zimmermann,
   LORIA/INRIA Nancy - Grand-Est.

   This file is part of the DPE Library.

   The DPE Library is free software; you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published by
   the Free Software Foundation; either version 2.1 of the License, or (at your
   option) any later version.

   The DPE Library is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
   License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with the DPE Library; see the file COPYING.LIB.  If not, write to
   the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
   MA 02110-1301, USA. */

/* WARNING: Patched version */

#ifndef FPLLL_DPE_H
#define FPLLL_DPE_H
#define __DPE

#include <math.h>   /* for round, floor, ceil */
#include <stdio.h>  /* For fprintf and fscanf */
#include <stdlib.h> /* For exit */
#if defined(__sun)  /* for round on Solaris 10 */
#include "ieeefp.h"
#include "tgmath.h"
#undef NAN
#define NAN 0.0 / 0.0
#endif
#include <limits.h>

#define DPE_VERSION_MAJOR 1
#define DPE_VERSION_MINOR 5

#if defined(__GNUC__) && (__GNUC__ >= 3)
#define DPE_LIKELY(x) (__builtin_expect(!!(x), 1))
#define DPE_UNLIKELY(x) (__builtin_expect((x), 0))
#define DPE_UNUSED_ATTR __attribute__((unused))
#else
#define DPE_LIKELY(x) (x)
#define DPE_UNLIKELY(x) (x)
#define DPE_UNUSED_ATTR
#endif

#if !defined(DPE_USE_DOUBLE) && !defined(DPE_USE_LONGDOUBLE)
#define DPE_USE_DOUBLE
#endif

#if (defined(__i386) || defined(__x86_64)) && !defined(DPE_LITTLEENDIAN32) &&                      \
    defined(DPE_USE_DOUBLE)
#define DPE_LITTLEENDIAN32
#endif

#if defined(DPE_USE_DOUBLE)
#define DPE_DOUBLE double /* mantissa type */
#define DPE_BITSIZE 53    /* bitsize of DPE_DOUBLE */
#define DPE_2_POW_BITSIZE 9007199254740992
/* DPE_LDEXP(DPE_DOUBLE m, DPEEXP e) return x = m * 2^e */
#define DPE_LDEXP ldexp
/* DPE_FREXP(DPE_DOUBLE x, DPEEXP *e) returns m, e such that x = m * 2^e with
   1/2 <= m < 1 */
#define DPE_FREXP frexp
/* DPE_ROUND(DPE_DOUBLE x) returns the nearest integer to x */
#define DPE_ROUND round
#define DPE_RINT rint
#define DPE_FLOOR floor
#define DPE_CEIL ceil
#define DPE_TRUNC trunc
#elif defined(DPE_USE_LONGDOUBLE)
#define DPE_DOUBLE long double
#define DPE_BITSIZE 64
#define DPE_2_POW_BITSIZE 18446744073709551616
#define DPE_LDEXP ldexpl
#define DPE_FREXP frexpl
#define DPE_ROUND roundl
#define DPE_RINT rintl
#define DPE_FLOOR floorl
#define DPE_CEIL ceill
#define DPE_TRUNC truncl
#else
#error "neither DPE_USE_DOUBLE nor DPE_USE_LONGDOUBLE is defined"
#endif

#if defined(DPE_USE_LONG)
#define DPE_EXP_T long      /* exponent type */
#define DPE_EXPMIN LONG_MIN /* smallest possible exponent */
#elif defined(DPE_USE_LONGLONG)
#define DPE_EXP_T long long
#define DPE_EXPMIN LONG_LONG_MIN
#else
#define DPE_EXP_T int      /* exponent type */
#define DPE_EXPMIN INT_MIN /* smallest possible exponent */
#endif

typedef union {
  double d;
  int i[2];
} dpe_double_words;

typedef struct
{
  DPE_DOUBLE d;  /* significand */
  DPE_EXP_T exp; /* exponent */
} dpe_struct;

typedef dpe_struct dpe_t[1];

#define DPE_MANT(x) ((x)->d)
#define DPE_EXP(x) ((x)->exp)
#define DPE_SIGN(x) ((DPE_MANT(x) < 0.0) ? -1 : (DPE_MANT(x) > 0.0))

#define DPE_INLINE static inline

/* initialize */
DPE_INLINE void dpe_init(dpe_t x DPE_UNUSED_ATTR) {}

/* clear */
DPE_INLINE void dpe_clear(dpe_t x DPE_UNUSED_ATTR) {}

/* set x to y */
DPE_INLINE void dpe_set(dpe_t x, const dpe_t y)
{
  DPE_MANT(x) = DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y);
}

/* set x to -y */
DPE_INLINE void dpe_neg(dpe_t x, const dpe_t y)
{
  DPE_MANT(x) = -DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y);
}

/* set x to |y| */
DPE_INLINE void dpe_abs(dpe_t x, const dpe_t y)
{
  DPE_MANT(x) = (DPE_MANT(y) >= 0) ? DPE_MANT(y) : -DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y);
}

/* set mantissa in [1/2, 1), except for 0 which has minimum exponent */
static void dpe_normalize(dpe_t x)
{
  if (DPE_UNLIKELY(DPE_MANT(x) == 0.0 || std::isfinite(DPE_MANT(x)) == 0))
  {
    if (DPE_MANT(x) == 0.0)
      DPE_EXP(x) = DPE_EXPMIN;
    /* otherwise let the exponent of NaN, Inf unchanged */
  }
  else
  {
    DPE_EXP_T e;
#ifdef DPE_LITTLEENDIAN32 /* 32-bit little endian */
    dpe_double_words dw;
    dw.d = DPE_MANT(x);
    e    = (dw.i[1] >> 20) & 0x7FF; /* unbiased exponent, 1022 for m=1/2 */
    DPE_EXP(x) += e - 1022;
    dw.i[1]     = (dw.i[1] & 0x800FFFFF) | 0x3FE00000;
    DPE_MANT(x) = dw.d;
#else /* portable code */
    double m    = DPE_MANT(x);
    DPE_MANT(x) = DPE_FREXP(m, &e);
    DPE_EXP(x) += e;
#endif
  }
}

#if defined(DPE_USE_DOUBLE)
static const double dpe_scale_tab[54] = {
    1.0000000000000000000,  0.5000000000000000000,  0.2500000000000000000,  0.1250000000000000000,
    0.0625000000000000000,  0.0312500000000000000,  0.0156250000000000000,  0.0078125000000000000,
    0.0039062500000000000,  0.0019531250000000000,  0.0009765625000000000,  0.0004882812500000000,
    0.0002441406250000000,  0.0001220703125000000,  0.0000610351562500000,  0.0000305175781250000,
    0.0000152587890625000,  7.6293945312500000e-6,  3.8146972656250000e-6,  1.9073486328125000e-6,
    9.5367431640625000e-7,  4.7683715820312500e-7,  2.3841857910156250e-7,  1.1920928955078125e-7,
    5.9604644775390625e-8,  2.9802322387695312e-8,  1.4901161193847656e-8,  7.4505805969238281e-9,
    3.7252902984619141e-9,  1.8626451492309570e-9,  9.3132257461547852e-10, 4.6566128730773926e-10,
    2.3283064365386963e-10, 1.1641532182693481e-10, 5.8207660913467407e-11, 2.9103830456733704e-11,
    1.4551915228366852e-11, 7.2759576141834259e-12, 3.6379788070917130e-12, 1.8189894035458565e-12,
    9.0949470177292824e-13, 4.5474735088646412e-13, 2.2737367544323206e-13, 1.1368683772161603e-13,
    5.6843418860808015e-14, 2.8421709430404007e-14, 1.4210854715202004e-14, 7.1054273576010019e-15,
    3.5527136788005009e-15, 1.7763568394002505e-15, 8.8817841970012523e-16, 4.4408920985006262e-16,
    2.2204460492503131e-16, 1.1102230246251565e-16};
#endif

DPE_INLINE double dpe_scale(double d, int s)
{
/* -DPE_BITSIZE < s <= 0 and 1/2 <= d < 1 */
#if defined(DPE_USE_DOUBLE)
  return d * dpe_scale_tab[-s];
#else /* portable code */
  return DPE_LDEXP(d, s);
#endif
}

/* set x to y */
DPE_INLINE void dpe_set_d(dpe_t x, double y)
{
  DPE_MANT(x) = (DPE_DOUBLE)y;
  DPE_EXP(x)  = 0;
  dpe_normalize(x);
}

/* set x to y */
DPE_INLINE void dpe_set_ld(dpe_t x, long double y)
{
  DPE_MANT(x) = (DPE_DOUBLE)y;
  DPE_EXP(x)  = 0;
  dpe_normalize(x);
}

/* set x to y */
DPE_INLINE void dpe_set_ui(dpe_t x, unsigned long y)
{
  DPE_MANT(x) = (DPE_DOUBLE)y;
  DPE_EXP(x)  = 0;
  dpe_normalize(x);
}

/* set x to y */
DPE_INLINE void dpe_set_si(dpe_t x, long y)
{
  DPE_MANT(x) = (DPE_DOUBLE)y;
  DPE_EXP(x)  = 0;
  dpe_normalize(x);
}

DPE_INLINE long dpe_get_si(const dpe_t x)
{
  DPE_DOUBLE d = DPE_LDEXP(DPE_MANT(x), DPE_EXP(x));
  return (long)d;
}

DPE_INLINE unsigned long dpe_get_ui(const dpe_t x)
{
  DPE_DOUBLE d = DPE_LDEXP(DPE_MANT(x), DPE_EXP(x));
  return (d < 0.0) ? 0 : (unsigned long)d;
}

DPE_INLINE double dpe_get_d(const dpe_t x) { return DPE_LDEXP(DPE_MANT(x), DPE_EXP(x)); }

DPE_INLINE long double dpe_get_ld(const dpe_t x) { return DPE_LDEXP(DPE_MANT(x), DPE_EXP(x)); }

#ifdef __GMP_H__
/* set x to y */
DPE_INLINE void dpe_set_z(dpe_t x, const mpz_t y)
{
  long e;
  DPE_MANT(x) = mpz_get_d_2exp(&e, y);
  DPE_EXP(x)  = (DPE_EXP_T)e;
}

/* set x to y, rounded to nearest */
DPE_INLINE void dpe_get_z(mpz_t x, const dpe_t y)
{
  DPE_EXP_T ey = DPE_EXP(y);
  if (ey >= DPE_BITSIZE) /* y is an integer */
  {
    DPE_DOUBLE d = DPE_MANT(y) * DPE_2_POW_BITSIZE; /* d is an integer */
    mpz_set_d(x, d);                                /* should be exact */
    mpz_mul_2exp(x, x, (unsigned long)ey - DPE_BITSIZE);
  }
  else /* DPE_EXP(y) < DPE_BITSIZE */
  {
    if (DPE_UNLIKELY(ey < 0)) /* |y| < 1/2 */
      mpz_set_ui(x, 0);
    else
    {
      DPE_DOUBLE d = DPE_LDEXP(DPE_MANT(y), ey);
      mpz_set_d(x, (double)DPE_ROUND(d));
    }
  }
}

/* return e and x such that y = x*2^e */
DPE_INLINE mp_exp_t dpe_get_z_exp(mpz_t x, const dpe_t y)
{
  mpz_set_d(x, DPE_MANT(y) * DPE_2_POW_BITSIZE);
  return DPE_EXP(y) - DPE_BITSIZE;
}

#endif  //#ifdef __GMP_H__

#ifdef __MPFR_H
/* set x to y where y is mpfr_t. */
DPE_INLINE void dpe_get_mpfr(mpfr_t y, const dpe_t x, const mp_rnd_t rnd)
{
  mpfr_set_d(y, DPE_MANT(x), rnd);
  mpfr_mul_2si(y, y, DPE_EXP(x), rnd);
}

/* set y to x */
DPE_INLINE void dpe_set_mpfr(dpe_t x, const mpfr_t y)
{
  long exp;
  DPE_MANT(x) = mpfr_get_d_2exp(&exp, y, GMP_RNDN);
  DPE_EXP(x)  = exp;
  dpe_normalize(x);
}

#endif  // __MPFR_H

/* x <- y + z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_add(dpe_t x, const dpe_t y, const dpe_t z)
{
  if (DPE_UNLIKELY(DPE_EXP(y) > DPE_EXP(z) + DPE_BITSIZE))
    /* |z| < 1/2*ulp(y), thus o(y+z) = y */
    dpe_set(x, y);
  else if (DPE_UNLIKELY(DPE_EXP(z) > DPE_EXP(y) + DPE_BITSIZE))
    dpe_set(x, z);
  else
  {
    DPE_EXP_T d = DPE_EXP(y) - DPE_EXP(z); /* |d| <= DPE_BITSIZE */

    if (d >= 0)
    {
      DPE_MANT(x) = DPE_MANT(y) + dpe_scale(DPE_MANT(z), -d);
      DPE_EXP(x)  = DPE_EXP(y);
    }
    else
    {
      DPE_MANT(x) = DPE_MANT(z) + dpe_scale(DPE_MANT(y), d);
      DPE_EXP(x)  = DPE_EXP(z);
    }
    dpe_normalize(x);
  }
}

/* x <- y - z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_sub(dpe_t x, const dpe_t y, const dpe_t z)
{
  if (DPE_UNLIKELY(DPE_EXP(y) > DPE_EXP(z) + DPE_BITSIZE))
    /* |z| < 1/2*ulp(y), thus o(y-z) = y */
    dpe_set(x, y);
  else if (DPE_UNLIKELY(DPE_EXP(z) > DPE_EXP(y) + DPE_BITSIZE))
    dpe_neg(x, z);
  else
  {
    DPE_EXP_T d = DPE_EXP(y) - DPE_EXP(z); /* |d| <= DPE_BITSIZE */

    if (d >= 0)
    {
      DPE_MANT(x) = DPE_MANT(y) - dpe_scale(DPE_MANT(z), -d);
      DPE_EXP(x)  = DPE_EXP(y);
    }
    else
    {
      DPE_MANT(x) = dpe_scale(DPE_MANT(y), d) - DPE_MANT(z);
      DPE_EXP(x)  = DPE_EXP(z);
    }
    dpe_normalize(x);
  }
}

/* x <- y * z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_mul(dpe_t x, const dpe_t y, const dpe_t z)
{
  DPE_MANT(x) = DPE_MANT(y) * DPE_MANT(z);
  DPE_EXP(x)  = DPE_EXP(y) + DPE_EXP(z);
  dpe_normalize(x);
}

/* x <- y * z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_mul_d(dpe_t x, const dpe_t y, const double z)
{
  dpe_t zt;
  dpe_set_d(zt, z);
  DPE_MANT(x) = DPE_MANT(y) * DPE_MANT(zt);
  DPE_EXP(x)  = DPE_EXP(y) + DPE_EXP(zt);
  dpe_normalize(x);
}

/* x <- x + y * z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_addmul(dpe_t x, const dpe_t y, const dpe_t z)
{
  dpe_t yz;
  dpe_mul(yz, y, z);
  dpe_add(x, x, yz);
}

/* x <- x - y * z, assuming y and z are normalized, returns x normalized */
DPE_INLINE void dpe_submul(dpe_t x, const dpe_t y, const dpe_t z)
{
  dpe_t yz;
  dpe_mul(yz, y, z);
  dpe_sub(x, x, yz);
}

/* x <- sqrt(y), assuming y is normalized, returns x normalized */
DPE_INLINE void dpe_sqrt(dpe_t x, const dpe_t y)
{
  DPE_EXP_T ey = DPE_EXP(y);
  if (ey % 2)
  {
    /* since 1/2 <= my < 1, 1/4 <= my/2 < 1 */
    DPE_MANT(x) = sqrt(0.5 * DPE_MANT(y));
    DPE_EXP(x)  = (ey + 1) / 2;
  }
  else
  {
    DPE_MANT(x) = sqrt(DPE_MANT(y));
    DPE_EXP(x)  = ey / 2;
  }
}

/* x <- y / z, assuming y and z are normalized, returns x normalized.
   Assumes z is not zero. */
DPE_INLINE void dpe_div(dpe_t x, const dpe_t y, const dpe_t z)
{
  DPE_MANT(x) = DPE_MANT(y) / DPE_MANT(z);
  DPE_EXP(x)  = DPE_EXP(y) - DPE_EXP(z);
  dpe_normalize(x);
}

/* x <- y * z, assuming y normalized, returns x normalized */
DPE_INLINE void dpe_mul_ui(dpe_t x, const dpe_t y, unsigned long z)
{
  DPE_MANT(x) = DPE_MANT(y) * (DPE_DOUBLE)z;
  DPE_EXP(x)  = DPE_EXP(y);
  dpe_normalize(x);
}

/* x <- y / z, assuming y normalized, z non-zero, returns x normalized */
DPE_INLINE void dpe_div_ui(dpe_t x, const dpe_t y, unsigned long z)
{
  DPE_MANT(x) = DPE_MANT(y) / (DPE_DOUBLE)z;
  DPE_EXP(x)  = DPE_EXP(y);
  dpe_normalize(x);
}

/* x <- y * 2^e */
DPE_INLINE void dpe_mul_2si(dpe_t x, const dpe_t y, long e)
{
  DPE_MANT(x) = DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y) + (DPE_EXP_T)e;
}

/* x <- y * 2^e */
DPE_INLINE void dpe_mul_2exp(dpe_t x, const dpe_t y, unsigned long e)
{
  DPE_MANT(x) = DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y) + (DPE_EXP_T)e;
}

/* x <- y / 2^e */
DPE_INLINE void dpe_div_2exp(dpe_t x, const dpe_t y, unsigned long e)
{
  DPE_MANT(x) = DPE_MANT(y);
  DPE_EXP(x)  = DPE_EXP(y) - (DPE_EXP_T)e;
}

/* return e and x such that y = x*2^e (equality is not guaranteed if the 'long'
   type has fewer bits than the significand in dpe_t) */
DPE_INLINE DPE_EXP_T dpe_get_si_exp(long *x, const dpe_t y)
{
  if (sizeof(long) == 4) /* 32-bit word: long has 31 bits */
  {
    *x = (long)(DPE_MANT(y) * 2147483648.0);
    return DPE_EXP(y) - 31;
  }
  else if (sizeof(long) == 8) /* 64-bit word: long has 63 bits */
  {
    *x = (long)(DPE_MANT(y) * 9223372036854775808.0);
    return DPE_EXP(y) - 63;
  }
  else
  {
    fprintf(stderr, "Error, neither 32-bit nor 64-bit word\n");
    abort();
  }
}

static DPE_UNUSED_ATTR int dpe_str_prec = 16;
static int dpe_out_str(FILE *s, int base, const dpe_t x) DPE_UNUSED_ATTR;

static int dpe_out_str(FILE *s, int base, const dpe_t x)
{
  DPE_DOUBLE d = DPE_MANT(x);
  DPE_EXP_T e2 = DPE_EXP(x);
  int e10      = 0;
  char sign    = ' ';
  if (DPE_UNLIKELY(base != 10))
  {
    fprintf(stderr, "Error in dpe_out_str, only base 10 allowed\n");
    abort();
  }
  if (d == 0.0)
#ifdef DPE_USE_DOUBLE
    return fprintf(s, "%1.*f", dpe_str_prec, d);
#else
    return fprintf(s, "foo\n %1.*Lf", dpe_str_prec, d);
#endif
  if (d < 0)
  {
    d    = -d;
    sign = '-';
  }
  if (e2 > 0)
  {
    while (e2 > 0)
    {
      e2--;
      d *= 2.0;
      if (d >= 10.0)
      {
        d /= 10.0;
        e10++;
      }
    }
  }
  else /* e2 <= 0 */
  {
    while (e2 < 0)
    {
      e2++;
      d /= 2.0;
      if (d < 1.0)
      {
        d *= 10.0;
        e10--;
      }
    }
  }
#ifdef DPE_USE_DOUBLE
  return fprintf(s, "%c%1.*f*10^%d", sign, dpe_str_prec, d, e10);
#else
  return fprintf(s, "%c%1.*Lf*10^%d", sign, dpe_str_prec, d, e10);
#endif
}

static size_t dpe_inp_str(dpe_t x, FILE *s, int base) DPE_UNUSED_ATTR;

static size_t dpe_inp_str(dpe_t x, FILE *s, int base)
{
  size_t res;
  DPE_DOUBLE d;
  if (DPE_UNLIKELY(base != 10))
  {
    fprintf(stderr, "Error in dpe_out_str, only base 10 allowed\n");
    abort();
  }
#ifdef DPE_USE_DOUBLE
  res = fscanf(s, "%lf", &d);
#else
  res = fscanf(s, "%Lf", &d);
#endif
  dpe_set_d(x, d);
  return res;
}

DPE_INLINE void dpe_dump(const dpe_t x)
{
  dpe_out_str(stdout, 10, x);
  putchar('\n');
}

DPE_INLINE int dpe_zero_p(const dpe_t x) { return DPE_MANT(x) == 0; }

/* return a positive value if x > y
   a negative value if x < y
   and 0 otherwise (x=y). */
DPE_INLINE int dpe_cmp(const dpe_t x, const dpe_t y)
{
  int sx = DPE_SIGN(x);
  int d  = sx - DPE_SIGN(y);

  if (d != 0)
    return d;
  else if (DPE_EXP(x) > DPE_EXP(y))
    return (sx > 0) ? 1 : -1;
  else if (DPE_EXP(y) > DPE_EXP(x))
    return (sx > 0) ? -1 : 1;
  else /* DPE_EXP(x) = DPE_EXP(y) */
    return (DPE_MANT(x) < DPE_MANT(y)) ? -1 : (DPE_MANT(x) > DPE_MANT(y));
}

DPE_INLINE int dpe_cmp_d(const dpe_t x, double d)
{
  dpe_t y;
  dpe_set_d(y, d);
  return dpe_cmp(x, y);
}

DPE_INLINE int dpe_cmp_ui(const dpe_t x, unsigned long d)
{
  dpe_t y;
  dpe_set_ui(y, d);
  return dpe_cmp(x, y);
}

DPE_INLINE int dpe_cmp_si(const dpe_t x, long d)
{
  dpe_t y;
  dpe_set_si(y, d);
  return dpe_cmp(x, y);
}

/* set x to integer nearest to y */
DPE_INLINE void dpe_round(dpe_t x, const dpe_t y)
{
  if (DPE_EXP(y) < 0) /* |y| < 1/2 */
    dpe_set_ui(x, 0);
  else if (DPE_EXP(y) >= DPE_BITSIZE) /* y is an integer */
    dpe_set(x, y);
  else
  {
    DPE_DOUBLE d;
    d = DPE_LDEXP(DPE_MANT(y), DPE_EXP(y));
    dpe_set_d(x, DPE_ROUND(d));
  }
}

/* set x to integer nearest to y */
DPE_INLINE void dpe_rint(dpe_t x, const dpe_t y)
{
  if (DPE_EXP(y) < 0) /* |y| < 1/2 */
    dpe_set_ui(x, 0);
  else if (DPE_EXP(y) >= DPE_BITSIZE) /* y is an integer */
    dpe_set(x, y);
  else
  {
    DPE_DOUBLE d;
    d = DPE_LDEXP(DPE_MANT(y), DPE_EXP(y));
    dpe_set_d(x, DPE_RINT(d));
  }
}

/* set x to the fractional part of y, defined as y - trunc(y), thus the
   fractional part has absolute value in [0, 1), and same sign as y */
DPE_INLINE void dpe_frac(dpe_t x, const dpe_t y)
{
  /* If |y| is smaller than 1, keep it */
  if (DPE_EXP(y) <= 0)
    dpe_set(x, y);
  else if (DPE_EXP(y) >= DPE_BITSIZE) /* y is an integer */
    dpe_set_ui(x, 0);
  else
  {
    DPE_DOUBLE d;
    d = DPE_LDEXP(DPE_MANT(y), DPE_EXP(y));
    dpe_set_d(x, d - DPE_TRUNC(d));
  }
}

/* set x to largest integer <= y */
DPE_INLINE void dpe_floor(dpe_t x, const dpe_t y)
{
  if (DPE_EXP(y) <= 0) /* |y| < 1 */
  {
    if (DPE_SIGN(y) >= 0) /* 0 <= y < 1 */
      dpe_set_ui(x, 0);
    else /* -1 < y < 0 */
      dpe_set_si(x, -1);
  }
  else if (DPE_EXP(y) >= DPE_BITSIZE) /* y is an integer */
    dpe_set(x, y);
  else
  {
    DPE_DOUBLE d;
    d = DPE_LDEXP(DPE_MANT(y), DPE_EXP(y));
    dpe_set_d(x, DPE_FLOOR(d));
  }
}

/* set x to smallest integer >= y */
DPE_INLINE void dpe_ceil(dpe_t x, const dpe_t y)
{
  if (DPE_EXP(y) <= 0) /* |y| < 1 */
  {
    if (DPE_SIGN(y) > 0) /* 0 < y < 1 */
      dpe_set_ui(x, 1);
    else /* -1 < y <= 0 */
      dpe_set_si(x, 0);
  }
  else if (DPE_EXP(y) >= DPE_BITSIZE) /* y is an integer */
    dpe_set(x, y);
  else
  {
    DPE_DOUBLE d;
    d = DPE_LDEXP(DPE_MANT(y), DPE_EXP(y));
    dpe_set_d(x, DPE_CEIL(d));
  }
}

DPE_INLINE void dpe_swap(dpe_t x, dpe_t y)
{
  DPE_EXP_T i  = DPE_EXP(x);
  DPE_DOUBLE d = DPE_MANT(x);
  DPE_EXP(x)   = DPE_EXP(y);
  DPE_MANT(x)  = DPE_MANT(y);
  DPE_EXP(y)   = i;
  DPE_MANT(y)  = d;
}

/* Ugly hacks: No correct rounding guaranteed */

DPE_INLINE void dpe_ugly_log(dpe_t x, const dpe_t y)
{
  dpe_set_d(x, ((double)DPE_EXP(y)) * M_LN2 + log(DPE_MANT(y)));
}

DPE_INLINE void dpe_ugly_exp(dpe_t x, const dpe_t y)
{
  // printf ("## exp is %ld\n", DPE_EXP(y));
  dpe_set_d(x, exp(((double)DPE_MANT(y)) * pow(2, ((double)DPE_EXP(y)))));
}

/* More hacks */
/* x = y^k */
DPE_INLINE void dpe_pow_si(dpe_t x, const dpe_t y, const unsigned int k)
{
  DPE_MANT(x) = pow(DPE_MANT(y), k);
  DPE_EXP(x)  = DPE_EXP(y) * k;
  dpe_normalize(x);
}

#ifdef __MPFR_H
/* x = e^y */
DPE_INLINE void dpe_exponential(dpe_t x, const dpe_t y)
{
  /* floor(log(DBL_MAX)) = 709 */
  if (dpe_cmp_ui(y, 709) <= 0 && dpe_cmp_si(y, -709) >= 0)
    dpe_ugly_exp(x, y);
  else
  {
    mpfr_t t, s;
    mpfr_init(t);
    mpfr_init(s);
    mpfr_set_d(t, 2.0, GMP_RNDN);
    mpfr_pow_si(t, t, DPE_EXP(y), GMP_RNDN);
    mpfr_set_d(s, DPE_MANT(y), GMP_RNDN);
    mpfr_mul(t, t, s, GMP_RNDN);
    mpfr_exp(t, t, GMP_RNDN);
    dpe_set_mpfr(x, t);
    mpfr_clear(t);
    mpfr_clear(s);
    mpfr_free_cache();
  }
}
#endif

DPE_INLINE void dpe_log(dpe_t x, const dpe_t y)
{
  dpe_set_d(x, (double)DPE_EXP(y) * log(2.0) + log(DPE_MANT(y)));
}

#endif /* FPLLL_DPE_H */