2
// INRIA Saclay - Ile de France (France).
3
// All rights reserved.
5
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
6
// modify it under the terms of the GNU Lesser General Public License as
7
// published by the Free Software Foundation; either version 3 of the License,
8
// or (at your option) any later version.
10
// Licensees holding a valid commercial license may use this file in
11
// accordance with the commercial license agreement provided with the software.
13
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
14
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16
// Author(s) : Marc Glisse
29
# include <CGAL/gmpxx.h>
31
# include <CGAL/gmp.h>
33
#include <CGAL/enum.h>
34
#include <CGAL/Interval_nt.h>
35
#include <CGAL/Gmpz.h>
36
#include <CGAL/Gmpq.h>
37
//#include <CGAL/Gmpzf.h>
39
#include <CGAL/Coercion_traits.h>
41
// The following is currently assumed in several places. I hope I am not
42
// making too many other assumptions.
43
// * limbs are 64 bits
44
// * if using gcc, sizeof(long long)==8
45
// * mpn_neg(_n) exists
47
// * not too fancy endianness
48
#if !defined(CGAL_DO_NOT_USE_MPZF) \
49
&& __GNU_MP_VERSION * 10 + __GNU_MP_VERSION_MINOR >= 43 \
50
&& GMP_NUMB_BITS == 64
51
#define CGAL_HAS_MPZF 1
53
// GMP-4.3.* has a different name for mpn_neg.
55
#define mpn_neg mpn_neg_n
57
// GMP before 5.0 doesn't provide mpn_copyi.
59
#define mpn_copyi(dst, src, siz) std::copy((src), (src)+(siz), (dst))
62
#include <boost/cstdint.hpp>
66
#pragma intrinsic(_BitScanForward64)
67
#pragma intrinsic(_BitScanReverse64)
70
#if defined(BOOST_MSVC)
71
# pragma warning(push)
72
# pragma warning(disable:4146 4244 4267 4800)
73
// warning on - applied on unsigned number
74
// conversion with loss of data
75
// conversion with loss of data
76
// int to bool performance
79
#if defined(__GNUC__) && defined(__GNUC_MINOR__) \
80
&& (__GNUC__ * 100 + __GNUC_MINOR__) >= 408 \
81
&& __cplusplus >= 201103L
82
#define CGAL_CAN_USE_CXX11_THREAD_LOCAL
86
#ifdef CGAL_MPZF_NO_USE_CACHE
87
# ifdef CGAL_MPZF_USE_CACHE
88
# undef CGAL_MPZF_USE_CACHE
91
# if !defined(CGAL_MPZF_USE_CACHE) \
92
&& defined(CGAL_HAS_THREADS) \
93
&& !defined(CGAL_I_PROMISE_I_WONT_USE_MANY_THREADS)
94
# define CGAL_MPZF_USE_CACHE 1
98
#define CGAL_MPZF_USE_CACHE 1
100
// On a dataset provided by Andreas, replacing Gmpq with this type in
101
// Epick reduced the running time of the construction of a Delaunay
102
// triangulation by a factor larger than 6
104
#if !defined(CGAL_HAS_THREADS)
105
#define CGAL_MPZF_THREAD_LOCAL
106
#define CGAL_MPZF_TLS
107
#elif defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
108
#define CGAL_MPZF_THREAD_LOCAL thread_local
109
#define CGAL_MPZF_TLS thread_local
110
#elif defined(_MSC_VER)
111
#define CGAL_MPZF_THREAD_LOCAL __declspec(thread)
112
#define CGAL_MPZF_TLS
114
#define CGAL_MPZF_THREAD_LOCAL __thread
115
#define CGAL_MPZF_TLS
116
// Too bad for the others
119
namespace Mpzf_impl {
120
// Warning: these pools aren't generic at all!
123
template <class T, class = void> struct pool1 {
124
static T pop() { T ret = data.back(); data.pop_back(); return ret; }
125
static void push(T t) { data.push_back(t); }
126
static bool empty() { return data.empty(); }
127
static const int extra = 0;
129
// thread_local would be fine, but not __declspec(thread) and for most
130
// compilers not __thread, since data is not POD.
131
static std::vector<T> data;
133
template <class T, class D> std::vector<T> pool1<T,D>::data;
135
// Use an intrusive single-linked list instead (allocate one more limb and use
136
// it to store the pointer to next), the difference isn't that noticable (still
137
// the list wins). Neither is thread-safe (both can be with threadlocal, and
138
// the list can be with an atomic compare-exchange (never tried)). With gcc,
139
// TLS has a large effect on classes with constructor/destructor, but is free
140
// for a simple pointer. The slowdown is because of PR55812.
142
// Leaks at thread destruction
143
template <class T, class = void> struct pool2 {
144
static T pop() { T ret = data(); data() = ptr(data()); return ret; }
145
static void push(T t) { ptr(t) = data(); data() = t; }
146
static bool empty() { return data() == 0; }
147
static const int extra = 1; // TODO: handle the case where a pointer is larger than a mp_limb_t
149
BOOST_STATIC_ASSERT(sizeof(T) >= sizeof(T*));
151
static CGAL_MPZF_TLS T data_ = 0;
154
static T& ptr(T t) { t -= extra+1; return *reinterpret_cast<T*>(t); }
157
#if defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
158
template <class T, class = void> struct pool3 {
159
static T pop() { T ret = data(); data() = ptr(data()); return ret; }
160
static void push(T t) { ptr(t) = data(); data() = t; }
161
static bool empty() { return data() == 0; }
162
static const int extra = 1; // TODO: handle the case where a pointer is larger than a mp_limb_t
164
BOOST_STATIC_ASSERT(sizeof(T) >= sizeof(T*));
168
// Deallocate everything. As an alternative, we could store it in a
169
// global location, for re-use by a later thread.
171
delete[] (pop() - (extra + 1));
175
static thread_local cleaner obj;
178
static T& ptr(T t) { t -= extra+1; return *reinterpret_cast<T*>(t); }
183
template <class T, class = void> struct no_pool {
184
static T pop() { throw "Shouldn't be here!"; }
185
static void push(T t) { delete [] (t - (extra+1)); }
186
static bool empty() { return true; }
187
static const int extra = 0;
190
// Only used with an argument known not to be 0.
191
inline int ctz (boost::uint64_t x) {
194
_BitScanForward64(&ret, x);
197
// Assume long long is 64 bits
198
return __builtin_ctzll (x);
201
inline int clz (boost::uint64_t x) {
204
_BitScanReverse64(&ret, x);
205
return 63 - (int)ret;
207
return __builtin_clzll (x);
211
// In C++11, std::fill_n returns a pointer to the end, but in C++03,
213
inline mp_limb_t* fill_n_ptr(mp_limb_t* p, int n, int c) {
214
#if __cplusplus >= 201103L
215
return std::fill_n (p, n, c);
217
mp_limb_t* q = p + n;
219
//std::fill_n (p, n, c);
220
//memset (p, sizeof(mp_limb_t)*n, c);
224
} // namespace Mpzf_impl
226
#undef CGAL_MPZF_THREAD_LOCAL
230
// * make data==0 a valid state for number 0. Incompatible with the cache. I
231
// tried, and it doesn't seem to help (may even hurt a bit).
234
#ifdef CGAL_MPZF_USE_CACHE
235
// More experiments to determine the best value would be good. It likely
236
// depends on the usage. Note that pool2 is fast enough that a conditional
237
// cache slows it down. A purely static cache (crash if it isn't large
238
// enough) still wins by about 11% on the Delaunay_3 construction, but is
239
// more complicated to handle.
240
// Evaluating a polynomial in double will never require more than roughly
241
// (2100*degree) bits, or (33*degree) mp_limb_t, which is very small. I
242
// checked by including an array of 150 limbs in every Mpzf (that's where
243
// the 11% number comes from).
244
// BONUS: doing that is thread-safe!
245
static const unsigned int cache_size = 8;
247
//#if !defined(CGAL_HAS_THREADS) || defined(CGAL_I_PROMISE_I_WONT_USE_MANY_THREADS)
248
// typedef Mpzf_impl::pool2<mp_limb_t*,Mpzf> pool;
249
//#elif defined(CGAL_CAN_USE_CXX11_THREAD_LOCAL)
250
// typedef Mpzf_impl::pool3<mp_limb_t*,Mpzf> pool;
252
typedef Mpzf_impl::no_pool<mp_limb_t*,Mpzf> pool;
255
mp_limb_t* data_; /* data_[0] is never 0 (except possibly for 0). */
256
inline mp_limb_t*& data() { return data_; };
257
inline mp_limb_t const* data() const { return data_; };
259
#ifdef CGAL_MPZF_USE_CACHE
260
mp_limb_t cache[cache_size + 1];
262
int size; /* Number of relevant limbs in data_. */
263
int exp; /* The number is data_ (an integer) * 2 ^ (64 * exp). */
268
void init(unsigned mini=2){
269
#ifdef CGAL_MPZF_USE_CACHE
270
if (mini <= cache_size) {
271
cache[0] = cache_size;
277
data() = pool::pop();
278
if(data()[-1] >= mini) return; // TODO: when mini==2, no need to check
279
delete[] (data() - (pool::extra+1)); // too small, useless
282
data() = (new mp_limb_t[mini+(pool::extra+1)]) + (pool::extra+1);
286
while(*--data()==0); // in case we skipped final zeroes
287
#ifdef CGAL_MPZF_USE_CACHE
288
if (data() == cache) return;
295
Mpzf(allocate,int i) { init(i); }
299
static void clear_pool () {
300
while (!pool::empty())
301
delete[] (pool::pop() - (pool::extra + 1));
307
Mpzf(): size(0), exp(0) {
310
Mpzf& operator=(Mpzf const& x){
311
unsigned asize=std::abs(x.size);
312
if(asize==0) { exp=0; size=0; return *this; }
313
if(this==&x) return *this;
314
while(*--data()==0); // factor that code somewhere?
316
#ifdef CGAL_MPZF_USE_CACHE
319
delete[] (data() - pool::extra);
324
mpn_copyi(data(),x.data(),asize);
328
int asize=std::abs(x.size);
332
if(size!=0) mpn_copyi(data(),x.data(),asize);
334
#if !defined(CGAL_CFG_NO_CPP0X_RVALUE_REFERENCE) \
335
&& !defined(CGAL_MPZF_USE_CACHE)
336
Mpzf(Mpzf&& x):data_(x.data()),size(x.size),exp(x.exp){
337
x.init(); // yes, that's a shame...
341
Mpzf& operator=(Mpzf&& x){
342
std::swap(size,x.size);
344
std::swap(data(),x.data());
347
friend Mpzf operator-(Mpzf&& x){
348
Mpzf ret = std::move(x);
349
ret.size = -ret.size;
353
Mpzf(int i) : exp(0) {
354
// assume that int is smaller than mp_limb_t
356
if (i == 0) { size = 0; }
357
else if (i > 0) { size = 1; data()[0] = i; }
358
else /* (i < 0) */ { size =-1; data()[0] = -(mp_limb_t)i; }
359
// cast to mp_limb_t because -INT_MIN is undefined
361
Mpzf(unsigned int i) : exp(0) {
362
// assume that int is smaller than mp_limb_t
364
if (i == 0) { size = 0; }
365
else /* (i > 0) */ { size = 1; data()[0] = i; }
367
Mpzf(long i) : exp(0) {
368
// assume that long is smaller than mp_limb_t
370
if (i == 0) { size = 0; }
371
else if (i > 0) { size = 1; data()[0] = i; }
372
else /* (i < 0) */ { size =-1; data()[0] = -(mp_limb_t)i; }
373
// cast to mp_limb_t because -LONG_MIN is undefined
375
Mpzf(unsigned long i) : exp(0) {
376
// assume that long is smaller than mp_limb_t
378
if (i == 0) { size = 0; }
379
else /* (i > 0) */ { size = 1; data()[0] = i; }
383
using boost::uint64_t;
385
#ifdef CGAL_LITTLE_ENDIAN
386
struct { uint64_t man:52; uint64_t exp:11; uint64_t sig:1; } s;
387
#else /* CGAL_BIG_ENDIAN */
389
struct { uint64_t sig:1; uint64_t exp:11; uint64_t man:52; } s;
395
uint64_t dexp = u.s.exp;
396
CGAL_assertion_msg(dexp != 2047, "Creating an Mpzf from infinity or NaN.");
398
if (d == 0) { size=0; exp=0; return; }
399
else { // denormal number
404
m = (1LL<<52) | u.s.man;
406
int e1 = (int)dexp+13;
407
// FIXME: make it more general! But not slower...
408
BOOST_STATIC_ASSERT(GMP_NUMB_BITS == 64);
411
// 52+1023+13==17*64 ?
413
// This seems very slightly faster
414
if(Mpzf_impl::ctz(m)+e2>=64){
415
data()[0] = m >> (64-e2);
420
if(e2>11){ // Wrong test for denormals
421
data()[1] = m >> (64-e2);
428
mp_limb_t d0 = (m << e2) & GMP_NUMB_MASK;
430
if (e2 != 0) // shifting by 64 is UB
431
d1 = m >> (GMP_NUMB_BITS - e2);
448
if(u.s.sig) size=-size;
449
//assert(to_double()==IA_force_to_double(d));
452
#ifdef CGAL_USE_GMPXX
453
Mpzf(mpz_class const&z){
454
init_from_mpz_t(z.get_mpz_t());
458
init_from_mpz_t(z.mpz());
460
void init_from_mpz_t(mpz_t const z){
461
exp=mpz_scan1(z,0)/GMP_NUMB_BITS;
462
size=mpz_size(z)-exp;
464
mpn_copyi(data(),z->_mp_d+exp,size);
468
// For debug purposes only
470
//std::cout << "size: " << size << std::endl;
471
if(size==0) { std::cout << "zero\n"; return; }
472
if(size<0) std::cout << "- ";
473
int asize = std::abs(size);
474
std::cout << std::hex;
475
while(--asize>=0) { std::cout << data()[asize] << ' '; }
476
std::cout << std::dec << "exp " << exp << ' ';
477
std::cout << std::dec << "size " << size << ' ';
478
asize = std::abs(size);
479
std::cout << "double: " << std::ldexp((double)data()[asize-1],64*(exp+asize-1))*((size<0)?-1:1) << '\n';
483
friend int Mpzf_abscmp(Mpzf const&a, Mpzf const&b){
484
int asize=std::abs(a.size);
485
int bsize=std::abs(b.size);
486
// size==0 should mean exp==-infinity, like with double.
487
// Since it doesn't, test for it explicitly.
488
if (bsize == 0) return asize;
489
if (asize == 0) return -1;
492
if(ah!=bh) return ah-bh;
493
int minsize=(std::min)(asize,bsize);
494
const mp_limb_t* adata=a.data()+(asize-1);
495
const mp_limb_t* bdata=b.data()+(bsize-1);
496
for(int i=0;i<minsize;++i,--adata,--bdata){
497
const mp_limb_t aa=*adata;
498
const mp_limb_t bb=*bdata;
499
if(aa!=bb) return (aa<bb)?-1:1;
501
return asize-bsize; // this assumes that we get rid of trailing zeros...
503
friend int Mpzf_cmp (Mpzf const&a, Mpzf const&b){
504
if ((a.size ^ b.size) < 0) return (a.size < 0) ? -1 : 1;
505
int res = Mpzf_abscmp(a, b);
506
return (a.size < 0) ? -res : res;
508
friend bool operator<(Mpzf const&a, Mpzf const&b){
509
if((a.size ^ b.size) < 0) return a.size < 0;
510
return ((a.size < 0) ? Mpzf_abscmp(b, a) : Mpzf_abscmp(a, b)) < 0;
512
friend bool operator>(Mpzf const&a, Mpzf const&b){
515
friend bool operator>=(Mpzf const&a, Mpzf const&b){
518
friend bool operator<=(Mpzf const&a, Mpzf const&b){
521
friend bool operator==(Mpzf const&a, Mpzf const&b){
522
if (a.exp != b.exp || a.size != b.size) return false;
523
if (a.size == 0) return true;
524
return mpn_cmp(a.data(), b.data(), std::abs(a.size)) == 0;
526
friend bool operator!=(Mpzf const&a, Mpzf const&b){
530
static Mpzf aors(Mpzf const&a, Mpzf const&b, int bsize){
533
int size=std::abs(a.size);
537
if(size!=0) mpn_copyi(res.data(),a.data(),size);
542
int size=std::abs(bsize);
546
mpn_copyi(res.data(),b.data(),size);
549
if((asize^bsize)>=0){
551
int absasize=std::abs(asize);
552
int absbsize=std::abs(bsize);
553
const mp_limb_t* adata=a.data();
554
const mp_limb_t* bdata=b.data();
557
if(aexp<bexp){ res.exp=a.exp; aexp=0; bexp=b.exp-a.exp; }
558
else { res.exp=b.exp; aexp=a.exp-b.exp; bexp=0; }
559
res.init((std::max)(absasize+aexp,absbsize+bexp)+1);
560
mp_limb_t* rdata=res.data();
562
// TODO: if aexp>0, swap a and b so we don't repeat the code.
564
if(absasize<=bexp){ // no overlap
565
mpn_copyi(rdata, adata, absasize);
567
rdata=Mpzf_impl::fill_n_ptr(rdata,bexp-absasize,0);
568
mpn_copyi(rdata, bdata, absbsize);
569
res.size=absbsize+bexp;
570
if(bsize<0) res.size=-res.size;
573
mpn_copyi(rdata, adata, bexp);
581
if(absbsize<=aexp){ // no overlap
582
mpn_copyi(rdata, bdata, absbsize);
584
rdata=Mpzf_impl::fill_n_ptr(rdata,aexp-absbsize,0);
585
mpn_copyi(rdata, adata, absasize);
586
res.size=absasize+aexp;
587
if(asize<0) res.size=-res.size;
590
mpn_copyi(rdata, bdata, aexp);
597
if(absasize>=absbsize){
598
mp_limb_t carry=mpn_add(rdata,adata,absasize,bdata,absbsize);
602
rdata[absasize]=carry;
605
mp_limb_t carry=mpn_add(rdata,bdata,absbsize,adata,absasize);
609
rdata[absbsize]=carry;
612
// unnecessary if a.exp != b.exp
613
while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
614
if(bsize<0) res.size=-res.size;
620
int cmp=Mpzf_abscmp(a,b);
621
if(cmp==0){ res.init(); res.size=0; res.exp=0; return res; }
622
if(cmp<0) { x=&b; y=&a; std::swap(xsize, ysize); }
624
int absxsize=std::abs(xsize);
625
int absysize=std::abs(ysize);
626
const mp_limb_t* xdata=x->data();
627
const mp_limb_t* ydata=y->data();
630
if(xexp<yexp){ res.exp=xexp; yexp-=xexp; xexp=0; }
631
else { res.exp=yexp; xexp-=yexp; yexp=0; }
632
res.init((std::max)(absxsize+xexp,absysize+yexp)+1);
633
mp_limb_t* rdata=res.data();
636
if(0<yexp){ // must have overlap since x is larger
637
mpn_copyi(rdata, xdata, yexp);
644
if(absysize<=xexp){ // no overlap
645
mpn_neg(rdata, ydata, absysize); // assert that it returns 1
647
rdata=Mpzf_impl::fill_n_ptr(rdata,xexp-absysize,-1);
648
mpn_sub_1(rdata, xdata, absxsize, 1);
649
res.size=absxsize+xexp;
650
if(res.data()[res.size-1]==0) --res.size;
651
if(xsize<0) res.size=-res.size;
654
mpn_neg(rdata, ydata, xexp); // assert that it returns 1
659
carry1=true; // assumes no trailing zeros
662
CGAL_assertion_code( mp_limb_t carry= )
663
mpn_sub(rdata, xdata, absxsize, ydata, absysize);
665
CGAL_assertion_code( carry+= )
666
mpn_sub_1(rdata, rdata, absxsize, 1);
667
CGAL_assertion(carry==0);
669
while(/*res.size>0&&*/res.data()[res.size-1]==0) --res.size;
670
while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
671
if(xsize<0) res.size=-res.size;
677
friend Mpzf operator+(Mpzf const&a, Mpzf const&b){
678
return aors(a,b,b.size);
681
friend Mpzf operator-(Mpzf const&a, Mpzf const&b){
682
return aors(a,b,-b.size);
685
friend Mpzf operator*(Mpzf const&a, Mpzf const&b){
686
int asize=std::abs(a.size);
687
int bsize=std::abs(b.size);
689
Mpzf res(allocate(),siz);
690
if(asize==0||bsize==0){res.exp=0;res.size=0;return res;}
694
high = mpn_mul(res.data(),a.data(),asize,b.data(),bsize);
696
high = mpn_mul(res.data(),b.data(),bsize,a.data(),asize);
698
if(res.data()[0]==0) { ++res.data(); ++res.exp; --siz; }
699
res.size=((a.size^b.size)>=0)?siz:-siz;
703
friend Mpzf Mpzf_square(Mpzf const&a){
704
int asize=std::abs(a.size);
706
Mpzf res(allocate(),siz);
708
if(asize==0){res.size=0;return res;}
709
mpn_sqr(res.data(),a.data(),asize);
710
mp_limb_t high = res.data()[siz-1];
712
if(res.data()[0]==0) { ++res.data(); ++res.exp; --siz; }
717
friend Mpzf operator/(Mpzf const&a, Mpzf const&b){
719
int asize=std::abs(a.size);
720
int bsize=std::abs(b.size);
721
int siz=asize+2-bsize;
722
Mpzf res(allocate(),asize+2);
723
if(bsize==0){throw std::range_error("Division by zero");}
724
if(asize==0){res.exp=0;res.size=0;return res;}
727
const mp_limb_t *adata = a.data();
728
const mp_limb_t *bdata = b.data();
729
mp_limb_t *qp = res.data();
730
mp_limb_t *rp = qp + siz;
731
if(Mpzf_impl::ctz(adata[0]) >= Mpzf_impl::ctz(bdata[0])){ // Easy case
733
mpn_tdiv_qr(qp, rp, 0, adata, asize, bdata, bsize);
735
for (int i=0; i<bsize; ++i)
736
if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
739
else if(adata[-1]==0){ // We are lucky
740
--adata; ++asize; --res.exp;
741
mpn_tdiv_qr(qp, rp, 0, adata, asize, bdata, bsize);
743
for (int i=0; i<bsize; ++i)
744
if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
749
Mpzf a2(allocate(),asize+1);
751
mpn_copyi(a2.data()+1,a.data(),asize);
752
// No need to complete a2, we just want the buffer.
753
//a2.size=(a.size<0)?(a.size-1):(a.size+1);
755
mpn_tdiv_qr(qp, rp, 0, a2.data(), asize+1, bdata, bsize);
757
for (int i=0; i<bsize; ++i)
758
if (rp[i] != 0) throw std::logic_error("non exact Mpzf division");
761
while(/*res.size>0&&*/res.data()[res.size-1]==0) --res.size;
762
//while(/*res.size>0&&*/res.data()[0]==0){--res.size;++res.data();++res.exp;}
763
if((a.size^b.size)<0) res.size=-res.size;
767
friend Mpzf Mpzf_gcd(Mpzf const&a, Mpzf const&b){
769
if (a.size == 0) return b;
770
if (b.size == 0) return a;
771
int asize=std::abs(a.size);
772
int bsize=std::abs(b.size);
773
int atz=Mpzf_impl::ctz(a.data()[0]);
774
int btz=Mpzf_impl::ctz(b.data()[0]);
775
int rtz=(std::min)(atz,btz);
776
Mpzf tmp(allocate(), asize);
777
Mpzf res(allocate(), bsize);
779
mpn_rshift(tmp.data(), a.data(), asize, atz);
780
if(tmp.data()[asize-1]==0) --asize;
782
else { mpn_copyi(tmp.data(), a.data(), asize); }
784
mpn_rshift(res.data(), b.data(), bsize, btz);
785
if(res.data()[bsize-1]==0) --bsize;
787
else { mpn_copyi(res.data(), b.data(), bsize); }
788
res.exp = 0; // Pick b.exp? or the average? 0 helps return 1 more often.
790
res.size = mpn_gcd(res.data(), res.data(), bsize, tmp.data(), asize);
792
res.size = mpn_gcd(res.data(), tmp.data(), asize, res.data(), bsize);
794
mp_limb_t c = mpn_lshift(res.data(), res.data(), res.size, rtz);
795
if(c) { res.data()[res.size]=c; ++res.size; }
800
friend bool Mpzf_is_square(Mpzf const&x){
801
if (x.size < 0) return false;
802
if (x.size == 0) return true;
803
// Assume that GMP_NUMB_BITS is even.
804
return mpn_perfect_square_p (x.data(), x.size);
807
friend Mpzf Mpzf_sqrt(Mpzf const&x){
809
if (x.size < 0) throw std::range_error("Sqrt of negative number");
810
if (x.size == 0) return 0;
811
if (x.exp % 2 == 0) {
812
Mpzf res(allocate(), (x.size + 1) / 2);
814
res.size = (x.size + 1) / 2;
815
CGAL_assertion_code(mp_size_t rem=)
816
mpn_sqrtrem(res.data(), 0, x.data(), x.size);
817
CGAL_assertion(rem==0);
820
else if (x.data()[-1] == 0) {
821
Mpzf res(allocate(), (x.size + 2) / 2);
822
res.exp = (x.exp - 1) / 2;
823
res.size = (x.size + 2) / 2;
824
CGAL_assertion_code(mp_size_t rem=)
825
mpn_sqrtrem(res.data(), 0, x.data()-1, x.size+1);
826
CGAL_assertion(rem==0);
830
Mpzf res(allocate(), (x.size + 2) / 2);
831
res.exp = (x.exp - 1) / 2;
832
res.size = (x.size + 2) / 2;
833
CGAL_assertion_code(mp_size_t rem=)
834
mpn_sqrtrem(res.data(), 0, x.data(), x.size);
835
CGAL_assertion(rem==0);
836
mpn_lshift(res.data(), res.data(), res.size, GMP_NUMB_BITS / 2);
841
friend Mpzf operator-(Mpzf const&x){
843
ret.size = -ret.size;
846
Mpzf& operator+=(Mpzf const&x){ *this=*this+x; return *this; }
847
Mpzf& operator-=(Mpzf const&x){ *this=*this-x; return *this; }
848
Mpzf& operator*=(Mpzf const&x){ *this=*this*x; return *this; }
850
bool is_canonical () const {
851
if (size == 0) return true;
852
if (data()[0] == 0) return false;
853
if (data()[std::abs(size)-1] == 0) return false;
857
bool is_zero () const {
861
bool is_one () const {
862
return exp==0 && size==1 && data()[0]==1;
865
CGAL::Sign sign () const { return CGAL::sign(size); }
867
double to_double () const {
868
// Assumes GMP_NUMB_BITS == 64
870
if(size==0) return 0;
871
int asize = std::abs(size);
872
mp_limb_t top = data()[asize-1];
874
if(top >= (1LL<<53) || asize == 1) /* ok */ ;
875
else { dtop += (double)data()[asize-2] * ldexp(1.,-GMP_NUMB_BITS); }
876
return ldexp( (size<0) ? -dtop : dtop, (asize-1+exp) * GMP_NUMB_BITS);
879
std::pair<double, double> to_interval () const {
880
// Assumes GMP_NUMB_BITS == 64
881
if (size == 0) return std::make_pair(0., 0.);
883
int asize = std::abs(size);
884
int e = 64 * (asize - 1 + exp);
885
mp_limb_t x = data()[asize-1];
886
int lz = Mpzf_impl::clz(x);
894
// Check for the few cases where dh=x works (asize==1 and the evicted
895
// bits from x were 0s)
897
else if (asize == 1) {
898
dl = dh = x; // conversion is exact
901
mp_limb_t y = data()[asize-2];
908
// Check for the few cases where dh=x works (asize==2 and the evicted
909
// bits from y were 0s)
911
typedef Interval_nt<> IA;
913
res = ldexp (res, e);
914
if (size < 0) res = -res;
915
return CGAL::to_interval(res);
916
// Use ldexp(Interval_nt,int) to delegate the hard thinking
917
// about over/underflow.
920
#ifdef CGAL_USE_GMPXX
921
#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
924
operator mpq_class () const {
926
export_to_mpq_t(q.get_mpq_t());
931
#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
934
operator Gmpq () const {
936
export_to_mpq_t(q.mpq());
939
void export_to_mpq_t(mpq_t q) const {
940
/* q must be 0/1 before this call */
941
CGAL_precondition(mpq_cmp_ui(q,0,1)==0);
943
mpz_import (mpq_numref (q),
945
-1, // data()[0] is the least significant part
947
0, // native endianness inside mp_limb_t
948
GMP_NAIL_BITS, // should be 0
951
mpq_mul_2exp(q, q, (sizeof(mp_limb_t) * CHAR_BIT * exp));
953
mpq_div_2exp(q, q, (sizeof(mp_limb_t) * CHAR_BIT * -exp));
960
#ifndef CGAL_CFG_NO_CPP0X_EXPLICIT_CONVERSION_OPERATORS
963
// This makes Mpzf==int ambiguous
964
operator Gmpzf () const {
966
z->_mp_d=const_cast<mp_limb_t*>(data());
969
// Only works for a very limited range of exponents
970
Gmpzf e(std::ldexp(1.,GMP_NUMB_BITS*exp));
975
friend void simplify_quotient(Mpzf& a, Mpzf& b){
976
// Avoid quotient(2^huge_a/2^huge_b)
979
// Simplify with gcd?
983
// Copied from Gmpzf, not sure that's the best thing to do.
985
std::ostream& operator<< (std::ostream& os, const Mpzf& a)
987
return os << a.to_double();
991
std::istream& operator>> (std::istream& is, Mpzf& a)
1001
template <> struct Algebraic_structure_traits< Mpzf >
1002
: public Algebraic_structure_traits_base< Mpzf, Integral_domain_without_division_tag > {
1003
typedef Tag_true Is_exact;
1004
typedef Tag_false Is_numerical_sensitive;
1007
: public std::unary_function< Type, bool > {
1008
bool operator()( const Type& x ) const {
1014
: public std::unary_function< Type, bool > {
1015
bool operator()( const Type& x ) const {
1021
: public std::binary_function< Type, Type, Type > {
1024
const Type& y ) const {
1025
return Mpzf_gcd(x, y);
1030
: public std::unary_function< Type, Type > {
1031
Type operator()( const Type& x ) const {
1032
return Mpzf_square(x);
1036
struct Integral_division
1037
: public std::binary_function< Type, Type, Type > {
1040
const Type& y ) const {
1046
: public std::unary_function< Type, Type > {
1047
Type operator()( const Type& x) const {
1048
return Mpzf_sqrt(x);
1053
: public std::binary_function< Type, Type&, bool > {
1054
bool operator()( const Type& x, Type& y ) const {
1055
// TODO: avoid doing 2 calls.
1056
if (!Mpzf_is_square(x)) return false;
1060
bool operator()( const Type& x) const {
1061
return Mpzf_is_square(x);
1066
template <> struct Real_embeddable_traits< Mpzf >
1067
: public INTERN_RET::Real_embeddable_traits_base< Mpzf , CGAL::Tag_true > {
1069
: public std::unary_function< Type, ::CGAL::Sign > {
1070
::CGAL::Sign operator()( const Type& x ) const {
1076
: public std::unary_function< Type, double > {
1077
double operator()( const Type& x ) const {
1078
return x.to_double();
1083
: public std::binary_function< Type, Type, Comparison_result > {
1084
Comparison_result operator()(
1086
const Type& y ) const {
1087
return CGAL::sign(Mpzf_cmp(x,y));
1092
: public std::unary_function< Type, std::pair< double, double > > {
1093
std::pair<double, double> operator()( const Type& x ) const {
1094
return x.to_interval();
1100
CGAL_DEFINE_COERCION_TRAITS_FOR_SELF(Mpzf)
1101
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short ,Mpzf)
1102
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int ,Mpzf)
1103
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long ,Mpzf)
1104
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float ,Mpzf)
1105
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(double ,Mpzf)
1106
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(Gmpz ,Mpzf)
1107
#ifdef CGAL_USE_GMPXX
1108
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(mpz_class,Mpzf)
1113
#if defined(BOOST_MSVC)
1114
# pragma warning(pop)
1117
#endif // GMP_NUMB_BITS == 64
1118
#endif // CGAL_MPZF_H
added
removed
include/CGAL/Mpzf.h
