libcgal-dev 4.9-1+b2 / usr / include / CGAL / Mpzf.h

// Copyright (c) 2013
// INRIA Saclay - Ile de France (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); 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 3 of the License,
// or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// Author(s)	:  Marc Glisse

#ifndef CGAL_MPZF_H
#define CGAL_MPZF_H
#include <cstdlib>
#include <algorithm>
#include <climits>
#include <vector>
#include <math.h>
#include <cmath>
#include <iostream>
#include <stdexcept>
#ifdef CGAL_USE_GMPXX
# include <CGAL/gmpxx.h>
#else
# include <CGAL/gmp.h>
#endif
#include <CGAL/enum.h>
#include <CGAL/Interval_nt.h>
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
//#include <CGAL/Gmpzf.h>

#include <CGAL/Coercion_traits.h>

// The following is currently assumed in several places. I hope I am not
// making too many other assumptions.
// * limbs are 64 bits
// * if using gcc, sizeof(long long)==8
// * mpn_neg(_n) exists
// * IEEE double
// * not too fancy endianness
#if !defined(CGAL_DO_NOT_USE_MPZF) \
    && __GNU_MP_VERSION * 10 + __GNU_MP_VERSION_MINOR >= 43 \
    && GMP_NUMB_BITS == 64
#define CGAL_HAS_MPZF 1

// GMP-4.3.* has a different name for mpn_neg.
#ifndef mpn_neg
#define mpn_neg mpn_neg_n
#endif
// GMP-4.3.0 is missing mpn_sqr.
#ifndef mpn_sqr
#define mpn_sqr(dest,a,n) mpn_mul_n(dest,a,a,n)
#endif
// GMP before 5.0 doesn't provide mpn_copyi.
#ifndef mpn_copyi
#define mpn_copyi(dst, src, siz) std::copy((src), (src)+(siz), (dst))
#endif

#include <boost/cstdint.hpp>

#ifdef _MSC_VER
#include <intrin.h>
#pragma intrinsic(_BitScanForward64)
#pragma intrinsic(_BitScanReverse64)
#endif

#ifdef __xlC__
#include <builtins.h>
#endif

#include <CGAL/assertions.h>
#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>
#include <boost/version.hpp>

#if defined(BOOST_MSVC)
#  pragma warning(push)
#  pragma warning(disable:4146 4244 4267 4800)
     // warning on - applied on unsigned number
     // conversion with loss of data
     // conversion with loss of data
     // int to bool performance
#endif


/*
#ifdef CGAL_MPZF_NO_USE_CACHE
# ifdef CGAL_MPZF_USE_CACHE
#  undef CGAL_MPZF_USE_CACHE
# endif
#else
# if !defined(CGAL_MPZF_USE_CACHE) \
     && defined(CGAL_HAS_THREADS) \
     && !defined(CGAL_I_PROMISE_I_WONT_USE_MANY_THREADS)
#  define CGAL_MPZF_USE_CACHE 1
# endif
#endif
*/
#define CGAL_MPZF_USE_CACHE 1

// On a dataset provided by Andreas, replacing Gmpq with this type in
// Epick reduced the running time of the construction of a Delaunay
// triangulation by a factor larger than 6

#if !defined(CGAL_HAS_THREADS)
#define CGAL_MPZF_THREAD_LOCAL
#define CGAL_MPZF_TLS
#elif defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
#define CGAL_MPZF_THREAD_LOCAL thread_local
#define CGAL_MPZF_TLS thread_local
#elif defined(_MSC_VER)
#define CGAL_MPZF_THREAD_LOCAL __declspec(thread)
#define CGAL_MPZF_TLS
#else
#define CGAL_MPZF_THREAD_LOCAL __thread
#define CGAL_MPZF_TLS
// Too bad for the others
#endif
namespace CGAL {
namespace Mpzf_impl {
// Warning: these pools aren't generic at all!

// Not thread-safe
template <class T, class = void> struct pool1 {
  static T pop() { T ret = data.back(); data.pop_back(); return ret; }
  static void push(T t) { data.push_back(t); }
  static bool empty() { return data.empty(); }
  static const int extra = 0;
  private:
  // thread_local would be fine, but not __declspec(thread) and for most
  // compilers not __thread, since data is not POD.
  static std::vector<T> data;
};
template <class T, class D> std::vector<T> pool1<T,D>::data;

// Use an intrusive single-linked list instead (allocate one more limb and use
// it to store the pointer to next), the difference isn't that noticable (still
// the list wins).  Neither is thread-safe (both can be with threadlocal, and
// the list can be with an atomic compare-exchange (never tried)).  With gcc,
// TLS has a large effect on classes with constructor/destructor, but is free
// for a simple pointer.  The slowdown is because of PR55812.

// Leaks at thread destruction
template <class T, class = void> struct pool2 {
  static T pop() { T ret = data(); data() = ptr(data()); return ret; }
  static void push(T t) { ptr(t) = data(); data() = t; }
  static bool empty() { return data() == 0; }
  static const int extra = 1; // TODO: handle the case where a pointer is larger than a mp_limb_t
  private:
  CGAL_static_assertion(sizeof(T) >= sizeof(T*));
  static T& data () {
    static CGAL_MPZF_TLS T data_ = 0;
    return data_;
  }
  static T& ptr(T t) { t -= extra+1; return *reinterpret_cast<T*>(t); }
};

#if defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
template <class T, class = void> struct pool3 {
  static T pop() { T ret = data(); data() = ptr(data()); return ret; }
  static void push(T t) { ptr(t) = data(); data() = t; }
  static bool empty() { return data() == 0; }
  static const int extra = 1; // TODO: handle the case where a pointer is larger than a mp_limb_t
  private:
  CGAL_static_assertion(sizeof(T) >= sizeof(T*));
  struct cleaner {
    T data_ = 0;
    ~cleaner(){
      // Deallocate everything. As an alternative, we could store it in a
      // global location, for re-use by a later thread.
      while (!empty())
	delete[] (pop() - (extra + 1));
    }
  };
  static T& data () {
    static thread_local cleaner obj;
    return obj.data_;
  }
  static T& ptr(T t) { t -= extra+1; return *reinterpret_cast<T*>(t); }
};
#endif

// No caching
template <class T, class = void> struct no_pool {
  static T pop() { throw "Shouldn't be here!"; }
  static void push(T t) { delete [] (t - (extra+1)); }
  static bool empty() { return true; }
  static const int extra = 0;
};

// Only used with an argument known not to be 0.
inline int ctz (boost::uint64_t x) {
#if defined(_MSC_VER)
  unsigned long ret;
  _BitScanForward64(&ret, x);
  return (int)ret;
#elif defined(__xlC__)
  return __cnttz8 (x);
#else
  // Assume long long is 64 bits
  return __builtin_ctzll (x);
#endif
}
inline int clz (boost::uint64_t x) {
#if defined(_MSC_VER)
  unsigned long ret;
  _BitScanReverse64(&ret, x);
  return 63 - (int)ret;
#elif defined(__xlC__)
  // Macro supposedly not defined on z/OS.
  return __cntlz8 (x);
#else
  return __builtin_clzll (x);
#endif
}

// In C++11, std::fill_n returns a pointer to the end, but in C++03,
// it returns void.
inline mp_limb_t* fill_n_ptr(mp_limb_t* p, int n, int c) {
#if __cplusplus >= 201103L
  return std::fill_n (p, n, c);
#else
  mp_limb_t* q = p + n;
  std::fill (p, q, c);
  //std::fill_n (p, n, c);
  //memset (p, sizeof(mp_limb_t)*n, c);
  return q;
#endif
}
} // namespace Mpzf_impl

#undef CGAL_MPZF_THREAD_LOCAL
#undef CGAL_MPZF_TLS

// TODO:
// * make data==0 a valid state for number 0. Incompatible with the cache. I
//   tried, and it doesn't seem to help (may even hurt a bit).
struct Mpzf {
  private:
#ifdef CGAL_MPZF_USE_CACHE
  // More experiments to determine the best value would be good. It likely
  // depends on the usage. Note that pool2 is fast enough that a conditional
  // cache slows it down. A purely static cache (crash if it isn't large
  // enough) still wins by about 11% on the Delaunay_3 construction, but is
  // more complicated to handle.
  // Evaluating a polynomial in double will never require more than roughly
  // (2100*degree) bits, or (33*degree) mp_limb_t, which is very small. I
  // checked by including an array of 150 limbs in every Mpzf (that's where
  // the 11% number comes from).
  // BONUS: doing that is thread-safe!
  static const unsigned int cache_size = 8;
#endif
//#if !defined(CGAL_HAS_THREADS) || defined(CGAL_I_PROMISE_I_WONT_USE_MANY_THREADS)
//  typedef Mpzf_impl::pool2<mp_limb_t*,Mpzf> pool;
//#elif defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
//  typedef Mpzf_impl::pool3<mp_limb_t*,Mpzf> pool;
//#else
  typedef Mpzf_impl::no_pool<mp_limb_t*,Mpzf> pool;
//#endif

  mp_limb_t* data_; /* data_[0] is never 0 (except possibly for 0). */
  inline mp_limb_t*& data() { return data_; };
  inline mp_limb_t const* data() const { return data_; };

#ifdef CGAL_MPZF_USE_CACHE
  mp_limb_t cache[cache_size + 1];
#endif
  int size; /* Number of relevant limbs in data_. */
  int exp; /* The number is data_ (an integer) * 2 ^ (64 * exp). */
  typedef int Exponent_type;
  typedef int Size_type;

  struct allocate{};
  struct noalloc{};

  void init(unsigned mini=2){
#ifdef CGAL_MPZF_USE_CACHE
    if (mini <= cache_size) {
      cache[0] = cache_size;
      data() = cache + 1;
      return;
    }
#endif
    if(!pool::empty()){
      data() = pool::pop();
      if(data()[-1] >= mini) return; // TODO: when mini==2, no need to check
      delete[] (data() - (pool::extra+1)); // too small, useless
    }
    if(mini<2) mini=2;
    data() = (new mp_limb_t[mini+(pool::extra+1)]) + (pool::extra+1);
    data()[-1] = mini;
  }
  void clear(){
    while(*--data()==0); // in case we skipped final zeroes
#ifdef CGAL_MPZF_USE_CACHE
    if (data() == cache) return;
#endif
    ++data();
    pool::push(data());
  }

  Mpzf(noalloc){}
  Mpzf(allocate,int i) { init(i); }

  public:

  static void clear_pool () {
    while (!pool::empty())
      delete[] (pool::pop() - (pool::extra + 1));
  }

  ~Mpzf(){
    clear();
  }
  Mpzf(): size(0), exp(0) {
    init();
  }
  Mpzf& operator=(Mpzf const& x){
    unsigned asize=std::abs(x.size);
    if(asize==0) { exp=0; size=0; return *this; }
    if(this==&x) return *this;
    while(*--data()==0); // factor that code somewhere?
    if(*data()<asize){
#ifdef CGAL_MPZF_USE_CACHE
      if (data() != cache)
#endif
	delete[] (data() - pool::extra);
      init(asize);
    } else ++data();
    size=x.size;
    exp=x.exp;
    mpn_copyi(data(),x.data(),asize);
    return *this;
  }
  Mpzf(Mpzf const& x){
    int asize=std::abs(x.size);
    init(asize);
    size=x.size;
    exp=x.exp;
    if(size!=0) mpn_copyi(data(),x.data(),asize);
  }
#if !defined(CGAL_CFG_NO_CPP0X_RVALUE_REFERENCE) \
    && !defined(CGAL_MPZF_USE_CACHE)
  Mpzf(Mpzf&& x):data_(x.data()),size(x.size),exp(x.exp){
    x.init(); // yes, that's a shame...
    x.size = 0;
    x.exp = 0;
  }
  Mpzf& operator=(Mpzf&& x){
    std::swap(size,x.size);
    exp = x.exp;
    std::swap(data(),x.data());
    return *this;
  }
  friend Mpzf operator-(Mpzf&& x){
    Mpzf ret = std::move(x);
    ret.size = -ret.size;
    return ret;
  }
#endif
  Mpzf(int i) : exp(0) {
    // assume that int is smaller than mp_limb_t
    init();
    if      (i == 0)    { size = 0; }
    else if (i >  0)    { size = 1; data()[0] = i; }
    else /* (i <  0) */ { size =-1; data()[0] = -(mp_limb_t)i; }
    // cast to mp_limb_t because -INT_MIN is undefined
  }
  Mpzf(unsigned int i) : exp(0) {
    // assume that int is smaller than mp_limb_t
    init();
    if      (i == 0)    { size = 0; }
    else /* (i >  0) */ { size = 1; data()[0] = i; }
  }
  Mpzf(long i) : exp(0) {
    // assume that long is smaller than mp_limb_t
    init();
    if      (i == 0)    { size = 0; }
    else if (i >  0)    { size = 1; data()[0] = i; }
    else /* (i <  0) */ { size =-1; data()[0] = -(mp_limb_t)i; }
    // cast to mp_limb_t because -LONG_MIN is undefined
  }
  Mpzf(unsigned long i) : exp(0) {
    // assume that long is smaller than mp_limb_t
    init();
    if      (i == 0)    { size = 0; }
    else /* (i >  0) */ { size = 1; data()[0] = i; }
  }
  Mpzf(double d){
    init();
    using boost::uint64_t;
    union {
#ifdef CGAL_LITTLE_ENDIAN
      struct { uint64_t man:52; uint64_t exp:11; uint64_t sig:1; } s;
#else /* CGAL_BIG_ENDIAN */
      //WARNING: untested!
      struct { uint64_t sig:1; uint64_t exp:11; uint64_t man:52; } s;
#endif
      double d;
    } u;
    u.d = d;
    uint64_t m;
    uint64_t dexp = u.s.exp;
    CGAL_assertion_msg(dexp != 2047, "Creating an Mpzf from infinity or NaN.");
    if (dexp == 0) {
      if (d == 0) { size=0; exp=0; return; }
      else { // denormal number
	m = u.s.man;
	++dexp;
      }
    } else {
      m = (1LL<<52) | u.s.man;
    }
    int e1 = (int)dexp+13;
    // FIXME: make it more general! But not slower...
    CGAL_static_assertion(GMP_NUMB_BITS == 64);
    int e2 = e1 % 64;
    exp = e1 / 64 - 17;
    // 52+1023+13==17*64 ?
#if 0
    // This seems very slightly faster
    if(Mpzf_impl::ctz(m)+e2>=64){
      data()[0] = m >> (64-e2);
      size = 1;
      ++exp;
    }else{
      data()[0] = m << e2;
      if(e2>11){ // Wrong test for denormals
	data()[1] = m >> (64-e2);
	size = 2;
      } else {
	size = 1;
      }
    }
#else
    mp_limb_t d0 = (m << e2) & GMP_NUMB_MASK;
    mp_limb_t d1 = 0;
    if (e2 != 0) // shifting by 64 is UB
      d1 = m >> (GMP_NUMB_BITS - e2);
    if (d0 == 0) {
      data()[0] = d1;
      size = 1;
      ++exp;
    }
    else {
      data()[0] = d0;
      if (d1 == 0) {
	size = 1;
      }
      else {
	data()[1] = d1;
	size = 2;
      }
    }
#endif
    if(u.s.sig) size=-size;
    //CGAL_assertion(to_double()==IA_force_to_double(d));
  }

#ifdef CGAL_USE_GMPXX
  Mpzf(mpz_class const&z){
    init_from_mpz_t(z.get_mpz_t());
  }
#endif
  Mpzf(Gmpz const&z){
    init_from_mpz_t(z.mpz());
  }
  void init_from_mpz_t(mpz_t const z){
    exp=Exponent_type(mpz_scan1(z,0)/GMP_NUMB_BITS);
    size=Size_type(mpz_size(z)-exp);
    init(size);
    mpn_copyi(data(),z->_mp_d+exp,size);
  }

#if 0
  // For debug purposes only
  void print()const{
    //std::cout << "size: " << size << std::endl;
    if(size==0) { std::cout << "zero\n"; return; }
    if(size<0) std::cout << "- ";
    int asize = std::abs(size);
    std::cout << std::hex;
    while(--asize>=0) { std::cout << data()[asize] << ' '; }
    std::cout << std::dec << "exp " << exp << ' ';
    std::cout << std::dec << "size " << size << ' ';
    asize = std::abs(size);
    std::cout << "double: " << std::ldexp((double)data()[asize-1],64*(exp+asize-1))*((size<0)?-1:1) << '\n';
  }
#endif

  friend int Mpzf_abscmp(Mpzf const&a, Mpzf const&b){
    int asize=std::abs(a.size);
    int bsize=std::abs(b.size);
    // size==0 should mean exp==-infinity, like with double.
    // Since it doesn't, test for it explicitly.
    if (bsize == 0) return asize;
    if (asize == 0) return -1;
    int ah=asize+a.exp;
    int bh=bsize+b.exp;
    if(ah!=bh) return ah-bh;
    int minsize=(std::min)(asize,bsize);
    const mp_limb_t* adata=a.data()+(asize-1);
    const mp_limb_t* bdata=b.data()+(bsize-1);
    for(int i=0;i<minsize;++i,--adata,--bdata){
      const mp_limb_t aa=*adata;
      const mp_limb_t bb=*bdata;
      if(aa!=bb) return (aa<bb)?-1:1;
    }
    return asize-bsize; // this assumes that we get rid of trailing zeros...
  }
  friend int Mpzf_cmp (Mpzf const&a, Mpzf const&b){
    if ((a.size ^ b.size) < 0) return (a.size < 0) ? -1 : 1;
    int res = Mpzf_abscmp(a, b);
    return (a.size < 0) ? -res : res;
  }
  friend bool operator<(Mpzf const&a, Mpzf const&b){
    if((a.size ^ b.size) < 0) return a.size < 0;
    return ((a.size < 0) ? Mpzf_abscmp(b, a) : Mpzf_abscmp(a, b)) < 0;
  }
  friend bool operator>(Mpzf const&a, Mpzf const&b){
    return b<a;
  }
  friend bool operator>=(Mpzf const&a, Mpzf const&b){
    return !(a<b);
  }
  friend bool operator<=(Mpzf const&a, Mpzf const&b){
    return !(a>b);
  }
  friend bool operator==(Mpzf const&a, Mpzf const&b){
    if (a.exp != b.exp || a.size != b.size) return false;
    if (a.size == 0) return true;
    return mpn_cmp(a.data(), b.data(), std::abs(a.size)) == 0;
  }
  friend bool operator!=(Mpzf const&a, Mpzf const&b){
    return !(a==b);
  }
  private:
  static Mpzf aors(Mpzf const&a, Mpzf const&b, int bsize){
    Mpzf res=noalloc();
    if(bsize==0){
      int size=std::abs(a.size);
      res.init(size);
      res.exp=a.exp;
      res.size=a.size;
      if(size!=0) mpn_copyi(res.data(),a.data(),size);
      return res;
    }
    int asize=a.size;
    if(asize==0){
      int size=std::abs(bsize);
      res.init(size);
      res.exp=b.exp;
      res.size=bsize;
      mpn_copyi(res.data(),b.data(),size);
      return res;
    }
    if((asize^bsize)>=0){
      // Addition
      int absasize=std::abs(asize);
      int absbsize=std::abs(bsize);
      const mp_limb_t* adata=a.data();
      const mp_limb_t* bdata=b.data();
      int aexp=a.exp;
      int bexp=b.exp;
      if(aexp<bexp){ res.exp=a.exp; aexp=0; bexp=b.exp-a.exp; }
      else { res.exp=b.exp; aexp=a.exp-b.exp; bexp=0; }
      res.init((std::max)(absasize+aexp,absbsize+bexp)+1);
      mp_limb_t* rdata=res.data();
      res.size=0;
      // TODO: if aexp>0, swap a and b so we don't repeat the code.
      if(0<bexp){
	if(absasize<=bexp){ // no overlap
	  mpn_copyi(rdata, adata, absasize);
	  rdata+=absasize;
	  rdata=Mpzf_impl::fill_n_ptr(rdata,bexp-absasize,0);
	  mpn_copyi(rdata, bdata, absbsize);
	  res.size=absbsize+bexp;
	  if(bsize<0) res.size=-res.size;
	  return res;
	} else {
	  mpn_copyi(rdata, adata, bexp);
	  adata+=bexp;
	  absasize-=bexp;
	  rdata+=bexp;
	  res.size=bexp;
	}
      }
      else if(0<aexp){
	if(absbsize<=aexp){ // no overlap
	  mpn_copyi(rdata, bdata, absbsize);
	  rdata+=absbsize;
	  rdata=Mpzf_impl::fill_n_ptr(rdata,aexp-absbsize,0);
	  mpn_copyi(rdata, adata, absasize);
	  res.size=absasize+aexp;
	  if(asize<0) res.size=-res.size;
	  return res;
	} else {
	  mpn_copyi(rdata, bdata, aexp);
	  bdata+=aexp;
	  absbsize-=aexp;
	  rdata+=aexp;
	  res.size=aexp;
	}
      }
      if(absasize>=absbsize){
	mp_limb_t carry=mpn_add(rdata,adata,absasize,bdata,absbsize);
	res.size+=absasize;
	if(carry!=0){
	  res.size++;
	  rdata[absasize]=carry;
	}
      } else {
	mp_limb_t carry=mpn_add(rdata,bdata,absbsize,adata,absasize);
	res.size+=absbsize;
	if(carry!=0){
	  res.size++;
	  rdata[absbsize]=carry;
	}
      }
      // unnecessary if a.exp != b.exp
      while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
      if(bsize<0) res.size=-res.size;
    } else {
      // Subtraction
      const Mpzf *x, *y;
      int xsize=a.size;
      int ysize=bsize;
      int cmp=Mpzf_abscmp(a,b);
      if(cmp==0){ res.init(); res.size=0; res.exp=0; return res; }
      if(cmp<0) { x=&b; y=&a; std::swap(xsize, ysize); }
      else { x=&a; y=&b; }
      int absxsize=std::abs(xsize);
      int absysize=std::abs(ysize);
      const mp_limb_t* xdata=x->data();
      const mp_limb_t* ydata=y->data();
      int xexp=x->exp;
      int yexp=y->exp;
      if(xexp<yexp){ res.exp=xexp; yexp-=xexp; xexp=0; }
      else { res.exp=yexp; xexp-=yexp; yexp=0; }
      res.init((std::max)(absxsize+xexp,absysize+yexp)+1);
      mp_limb_t* rdata=res.data();
      res.size=0;
      bool carry1=false;
      if(0<yexp){ // must have overlap since x is larger
	mpn_copyi(rdata, xdata, yexp);
	xdata+=yexp;
	absxsize-=yexp;
	rdata+=yexp;
	res.size=yexp;
      }
      else if(0<xexp){
	if(absysize<=xexp){ // no overlap
	  mpn_neg(rdata, ydata, absysize); // assert that it returns 1
	  rdata+=absysize;
	  rdata=Mpzf_impl::fill_n_ptr(rdata,xexp-absysize,-1);
	  mpn_sub_1(rdata, xdata, absxsize, 1);
	  res.size=absxsize+xexp;
	  if(res.data()[res.size-1]==0) --res.size;
	  if(xsize<0) res.size=-res.size;
	  return res;
	} else {
	  mpn_neg(rdata, ydata, xexp); // assert that it returns 1
	  ydata+=xexp;
	  absysize-=xexp;
	  rdata+=xexp;
	  res.size=xexp;
	  carry1=true; // assumes no trailing zeros
	}
      }
      CGAL_assertion_code( mp_limb_t carry= )
	mpn_sub(rdata, xdata, absxsize, ydata, absysize);
      if(carry1)
	CGAL_assertion_code( carry+= )
	  mpn_sub_1(rdata, rdata, absxsize, 1);
      CGAL_assertion(carry==0);
      res.size+=absxsize;
      while(/*res.size>0&&*/res.data()[res.size-1]==0) --res.size;
      while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
      if(xsize<0) res.size=-res.size;
    }
    return res;
  }

  public:
  friend Mpzf operator+(Mpzf const&a, Mpzf const&b){
    return aors(a,b,b.size);
  }

  friend Mpzf operator-(Mpzf const&a, Mpzf const&b){
    return aors(a,b,-b.size);
  }

  friend Mpzf operator*(Mpzf const&a, Mpzf const&b){
    int asize=std::abs(a.size);
    int bsize=std::abs(b.size);
    int siz=asize+bsize;
    Mpzf res(allocate(),siz);
    if(asize==0||bsize==0){res.exp=0;res.size=0;return res;}
    res.exp=a.exp+b.exp;
    mp_limb_t high;
    if(asize>=bsize)
      high = mpn_mul(res.data(),a.data(),asize,b.data(),bsize);
    else
      high = mpn_mul(res.data(),b.data(),bsize,a.data(),asize);
    if(high==0) --siz;
    if(res.data()[0]==0) { ++res.data(); ++res.exp; --siz; }
    res.size=((a.size^b.size)>=0)?siz:-siz;
    return res;
  }

  friend Mpzf Mpzf_square(Mpzf const&a){
    int asize=std::abs(a.size);
    int siz=2*asize;
    Mpzf res(allocate(),siz);
    res.exp=2*a.exp;
    if(asize==0){res.size=0;return res;}
    mpn_sqr(res.data(),a.data(),asize);
    mp_limb_t high = res.data()[siz-1];
    if(high==0) --siz;
    if(res.data()[0]==0) { ++res.data(); ++res.exp; --siz; }
    res.size=siz;
    return res;
  }

  friend Mpzf operator/(Mpzf const&a, Mpzf const&b){
    // FIXME: Untested
    int asize=std::abs(a.size);
    int bsize=std::abs(b.size);
    int siz=asize+2-bsize;
    Mpzf res(allocate(),asize+2);
    if(bsize==0){throw std::range_error("Division by zero");}
    if(asize==0){res.exp=0;res.size=0;return res;}
    res.size=siz;
    res.exp=a.exp-b.exp;
    const mp_limb_t *adata = a.data();
    const mp_limb_t *bdata = b.data();
    mp_limb_t *qp = res.data();
    mp_limb_t *rp = qp + siz;
    if(Mpzf_impl::ctz(adata[0]) >= Mpzf_impl::ctz(bdata[0])){ // Easy case
      --res.size;
      mpn_tdiv_qr(qp, rp, 0, adata, asize, bdata, bsize);
      CGAL_assertion_code(
	  for (int i=0; i<bsize; ++i)
	    if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
      )
    }
    else if(adata[-1]==0){ // We are lucky
      --adata; ++asize; --res.exp;
      mpn_tdiv_qr(qp, rp, 0, adata, asize, bdata, bsize);
      CGAL_assertion_code(
	  for (int i=0; i<bsize; ++i)
	    if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
      )
    }
    else{
      --res.exp;
      Mpzf a2(allocate(),asize+1);
      a2.data()[0]=0;
      mpn_copyi(a2.data()+1,a.data(),asize);
      // No need to complete a2, we just want the buffer.
      //a2.size=(a.size<0)?(a.size-1):(a.size+1);
      //a2.exp = a.exp-1;
      mpn_tdiv_qr(qp, rp, 0, a2.data(), asize+1, bdata, bsize);
      CGAL_assertion_code(
	  for (int i=0; i<bsize; ++i)
	    if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
      )
    }
    while(/*res.size>0&&*/res.data()[res.size-1]==0) --res.size;
    //while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
    if((a.size^b.size)<0) res.size=-res.size;
    return res;
  }

  friend Mpzf Mpzf_gcd(Mpzf const&a, Mpzf const&b){
    // FIXME: Untested
    if (a.size == 0) return b;
    if (b.size == 0) return a;
    int asize=std::abs(a.size);
    int bsize=std::abs(b.size);
    int atz=Mpzf_impl::ctz(a.data()[0]);
    int btz=Mpzf_impl::ctz(b.data()[0]);
    int rtz=(std::min)(atz,btz);
    Mpzf tmp(allocate(), asize);
    Mpzf res(allocate(), bsize);
    if (atz != 0) {
      mpn_rshift(tmp.data(), a.data(), asize, atz);
      if(tmp.data()[asize-1]==0) --asize;
    }
    else { mpn_copyi(tmp.data(), a.data(), asize); }
    if (btz != 0) {
      mpn_rshift(res.data(), b.data(), bsize, btz);
      if(res.data()[bsize-1]==0) --bsize;
    }
    else { mpn_copyi(res.data(), b.data(), bsize); }
    res.exp = 0; // Pick b.exp? or the average? 0 helps return 1 more often.
    if (asize < bsize)
      res.size = Size_type(mpn_gcd(res.data(), res.data(), bsize, tmp.data(), asize));
    else
      res.size = Size_type(mpn_gcd(res.data(), tmp.data(), asize, res.data(), bsize));
    if(rtz!=0) {
      mp_limb_t c = mpn_lshift(res.data(), res.data(), res.size, rtz);
      if(c) { res.data()[res.size]=c; ++res.size; }
    }
    return res;
  }

  friend bool Mpzf_is_square(Mpzf const&x){
    if (x.size < 0) return false;
    if (x.size == 0) return true;
    // Assume that GMP_NUMB_BITS is even.
    return mpn_perfect_square_p (x.data(), x.size);
  }

  friend Mpzf Mpzf_sqrt(Mpzf const&x){
    // FIXME: Untested
    if (x.size < 0) throw std::range_error("Sqrt of negative number");
    if (x.size == 0) return 0;
    if (x.exp % 2 == 0) {
      Mpzf res(allocate(), (x.size + 1) / 2);
      res.exp = x.exp / 2;
      res.size = (x.size + 1) / 2;
      CGAL_assertion_code(mp_size_t rem=)
      mpn_sqrtrem(res.data(), 0, x.data(), x.size);
      CGAL_assertion(rem==0);
      return res;
    }
    else if (x.data()[-1] == 0) {
      Mpzf res(allocate(), (x.size + 2) / 2);
      res.exp = (x.exp - 1) / 2;
      res.size = (x.size + 2) / 2;
      CGAL_assertion_code(mp_size_t rem=)
      mpn_sqrtrem(res.data(), 0, x.data()-1, x.size+1);
      CGAL_assertion(rem==0);
      return res;
    }
    else {
      Mpzf res(allocate(), (x.size + 2) / 2);
      res.exp = (x.exp - 1) / 2;
      res.size = (x.size + 2) / 2;
      CGAL_assertion_code(mp_size_t rem=)
      mpn_sqrtrem(res.data(), 0, x.data(), x.size);
      CGAL_assertion(rem==0);
      mpn_lshift(res.data(), res.data(), res.size, GMP_NUMB_BITS / 2);
      return res;
    }
  }

  friend Mpzf operator-(Mpzf const&x){
    Mpzf ret = x;
    ret.size = -ret.size;
    return ret;
  }
  Mpzf& operator+=(Mpzf const&x){ *this=*this+x; return *this; }
  Mpzf& operator-=(Mpzf const&x){ *this=*this-x; return *this; }
  Mpzf& operator*=(Mpzf const&x){ *this=*this*x; return *this; }
  Mpzf& operator/=(Mpzf const&x){ *this=*this/x; return *this; }

  bool is_canonical () const {
    if (size == 0) return true;
    if (data()[0] == 0) return false;
    if (data()[std::abs(size)-1] == 0) return false;
    return true;
  }

  bool is_zero () const {
    return size==0;
  }

  bool is_one () const {
    return exp==0 && size==1 && data()[0]==1;
  }

  CGAL::Sign sign () const { return CGAL::sign(size); }

  double to_double () const {
    // Assumes GMP_NUMB_BITS == 64
    using std::ldexp;
    if(size==0) return 0;
    int asize = std::abs(size);
    mp_limb_t top = data()[asize-1];
    double dtop = (double)top;
    if(top >= (1LL<<53) || asize == 1) /* ok */ ;
    else { dtop += (double)data()[asize-2] * ldexp(1.,-GMP_NUMB_BITS); }
    return ldexp( (size<0) ? -dtop : dtop, (asize-1+exp) * GMP_NUMB_BITS);
  }

  std::pair<double, double> to_interval () const {
    // Assumes GMP_NUMB_BITS == 64
    if (size == 0) return std::make_pair(0., 0.);
    double dl, dh;
    int asize = std::abs(size);
    int e = 64 * (asize - 1 + exp);
    mp_limb_t x = data()[asize-1];
    int lz = Mpzf_impl::clz(x);
    if (lz <= 11) {
      if (lz != 11) {
	e += (11 - lz);
	x >>= (11 - lz);
      }
      dl = double(x);
      dh = double(x + 1);
      // Check for the few cases where dh=x works (asize==1 and the evicted
      // bits from x were 0s)
    }
    else if (asize == 1) {
      dl = dh = double(x); // conversion is exact
    }
    else {
      mp_limb_t y = data()[asize-2];
      e -= (lz - 11);
      x <<= (lz - 11);
      y >>= (75 - lz);
      x |= y;
      dl = double(x);
      dh = double(x + 1);
      // Check for the few cases where dh=x works (asize==2 and the evicted
      // bits from y were 0s)
    }
    typedef Interval_nt<> IA;
    IA res (dl, dh);
    res = ldexp (res, e);
    if (size < 0) res = -res;
    return CGAL::to_interval(res);
    // Use ldexp(Interval_nt,int) to delegate the hard thinking
    // about over/underflow.
  }

#ifdef CGAL_USE_GMPXX
#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
  explicit
#endif
  operator mpq_class () const {
    mpq_class q;
    export_to_mpq_t(q.get_mpq_t());
    return q;
  }
#endif

#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
  explicit
#endif
  operator Gmpq () const {
    Gmpq q;
    export_to_mpq_t(q.mpq());
    return q;
  }
  void export_to_mpq_t(mpq_t q) const {
    /* q must be 0/1 before this call */
    CGAL_precondition(mpq_cmp_ui(q,0,1)==0);
    if (size != 0) {
      mpz_import (mpq_numref (q),
		  std::abs(size),
		  -1, // data()[0] is the least significant part
		  sizeof(mp_limb_t),
		  0, // native endianness inside mp_limb_t
		  GMP_NAIL_BITS, // should be 0
		  data());
      if (exp > 0)
	mpq_mul_2exp(q, q, (sizeof(mp_limb_t) * CHAR_BIT *  exp));
      else if (exp < 0)
	mpq_div_2exp(q, q, (sizeof(mp_limb_t) * CHAR_BIT * -exp));

      if (size < 0)
	mpq_neg(q,q);
    }
  }
#if 0
#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
  explicit
#endif
// This makes Mpzf==int ambiguous
  operator Gmpzf () const {
    mpz_t z;
    z->_mp_d=const_cast<mp_limb_t*>(data());
    z->_mp_size=size;
    Gmpzf m(z);
// Only works for a very limited range of exponents
    Gmpzf e(std::ldexp(1.,GMP_NUMB_BITS*exp));
    return m*e;
  }
#endif

  friend void simplify_quotient(Mpzf& a, Mpzf& b){
    // Avoid quotient(2^huge_a/2^huge_b)
    a.exp -= b.exp;
    b.exp = 0;
    // Simplify with gcd?
  }
};

// Copied from Gmpzf, not sure that's the best thing to do.
inline
std::ostream& operator<< (std::ostream& os, const Mpzf& a)
{
    return os << a.to_double();
}

inline
std::istream& operator>> (std::istream& is, Mpzf& a)
{
  double d;
  is >> d;
  if (is)
    a = d;
  return is;
}


  template <> struct Algebraic_structure_traits< Mpzf >
    : public Algebraic_structure_traits_base< Mpzf, Integral_domain_without_division_tag >  {
      typedef Tag_true            Is_exact;
      typedef Tag_false            Is_numerical_sensitive;

      struct Is_zero
	: public std::unary_function< Type, bool > {
	  bool operator()( const Type& x ) const {
	    return x.is_zero();
	  }
	};

      struct Is_one
	: public std::unary_function< Type, bool > {
	  bool operator()( const Type& x ) const {
	    return x.is_one();
	  }
	};

      struct Gcd
	: public std::binary_function< Type, Type, Type > {
	  Type operator()(
	      const Type& x,
	      const Type& y ) const {
	    return Mpzf_gcd(x, y);
	  }
	};

      struct Square
	: public std::unary_function< Type, Type > {
	  Type operator()( const Type& x ) const {
	    return Mpzf_square(x);
	  }
	};

      struct Integral_division
	: public std::binary_function< Type, Type, Type > {
	  Type operator()(
	      const Type& x,
	      const Type& y ) const {
	    return x / y;
	  }
	};

      struct Sqrt
	: public std::unary_function< Type, Type > {
	  Type operator()( const Type& x) const {
	    return Mpzf_sqrt(x);
	  }
	};

      struct Is_square
	: public std::binary_function< Type, Type&, bool > {
	  bool operator()( const Type& x, Type& y ) const {
	    // TODO: avoid doing 2 calls.
	    if (!Mpzf_is_square(x)) return false;
	    y = Mpzf_sqrt(x);
	    return true;
	  }
	  bool operator()( const Type& x) const {
	    return Mpzf_is_square(x);
	  }
	};

    };
  template <> struct Real_embeddable_traits< Mpzf >
    : public INTERN_RET::Real_embeddable_traits_base< Mpzf , CGAL::Tag_true > {
      struct Sgn
	: public std::unary_function< Type, ::CGAL::Sign > {
	  ::CGAL::Sign operator()( const Type& x ) const {
	    return x.sign();
	  }
	};

      struct To_double
	: public std::unary_function< Type, double > {
	    double operator()( const Type& x ) const {
	      return x.to_double();
	    }
	};

      struct Compare
	: public std::binary_function< Type, Type, Comparison_result > {
	    Comparison_result operator()(
		const Type& x,
		const Type& y ) const {
	      return CGAL::sign(Mpzf_cmp(x,y));
	    }
	};

      struct To_interval
	: public std::unary_function< Type, std::pair< double, double > > {
	    std::pair<double, double> operator()( const Type& x ) const {
	      return x.to_interval();
	    }
	};

    };

CGAL_DEFINE_COERCION_TRAITS_FOR_SELF(Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short    ,Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int      ,Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long     ,Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float    ,Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(double   ,Mpzf)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(Gmpz     ,Mpzf)
#ifdef CGAL_USE_GMPXX
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(mpz_class,Mpzf)
#endif

}

/* There isn't much Eigen can do with such a type,
 * mostly this is here for IsInteger to protect people.
 */
namespace Eigen {
  template<class> struct NumTraits;
  template<> struct NumTraits<CGAL::Mpzf>
  {
    typedef CGAL::Mpzf Real;
    /* Should this be Quotient<Mpzf>? Gmpq?  */
    typedef CGAL::Mpzf NonInteger;
    typedef CGAL::Mpzf Nested;
    typedef CGAL::Mpzf Literal;

    static inline Real epsilon() { return 0; }
    static inline Real dummy_precision() { return 0; }

    enum {
      /* Only exact divisions are supported, close enough to an integer.
       * This way we get compilation failures instead of runtime.  */
      IsInteger = 1,
      IsSigned = 1,
      IsComplex = 0,
      RequireInitialization = 1,
      ReadCost = 6,
      AddCost = 30,
      MulCost = 50
    };
  };
}

#if defined(BOOST_MSVC)
#  pragma warning(pop)
#endif

#endif // GMP_NUMB_BITS == 64
#endif // CGAL_MPZF_H