1 // Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
2 //
3 // This file is part of CGAL (www.cgal.org)
4 //
5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Algebraic_kernel_d/include/CGAL/RS/dyadic.h $
6 // $Id: dyadic.h 52164b1 2019-10-19T15:34:59+02:00 Sébastien Loriot
7 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
8 //
9 // Author: Luis Peñaranda <luis.penaranda@gmx.com>
10 
11 #ifndef CGAL_RS_DYADIC_H
12 #define CGAL_RS_DYADIC_H
13 
14 #include <stdio.h>
15 #include <math.h>
16 #include <gmp.h>
17 #include <mpfr.h>
18 #include <CGAL/assertions.h>
19 
20 // for c++, compile with -lgmpxx
21 #ifdef __cplusplus
22 #include <iostream>
23 #endif
24 
25 #define CGALRS_dyadic_struct            __mpfr_struct
26 #define CGALRS_dyadic_t                 mpfr_t
27 #define CGALRS_dyadic_ptr               mpfr_ptr
28 #define CGALRS_dyadic_srcptr            mpfr_srcptr
29 
30 // some auxiliary defines
31 #define CGALRS_dyadic_set_prec(D,P) \
32  ( mpfr_set_prec( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN) )
33 
34 #define CGALRS_dyadic_prec_round(D,P) \
35  ( mpfr_prec_round( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN, GMP_RNDN) )
36 
37 #define CGALRS_dyadic_set_exp(D,E) \
38  ( CGAL_assertion( (E) <= mpfr_get_emax() && \
39                    (E) >= mpfr_get_emin() ) ,\
40    mpfr_set_exp(D,E) )
41 
42 // init functions
43 #define CGALRS_dyadic_init(D)          mpfr_init2(D,MPFR_PREC_MIN)
44 #define CGALRS_dyadic_init2(D,P)       mpfr_init2(D,P)
45 #define CGALRS_dyadic_clear(D)         mpfr_clear(D)
46 
CGALRS_dyadic_set(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op)47 inline void CGALRS_dyadic_set(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
48         if(rop!=op){
49                 CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
50                 mpfr_set(rop,op,GMP_RNDN);
51         }
52         CGAL_assertion(mpfr_equal_p(rop,op)!=0);
53 }
54 
CGALRS_dyadic_set_z(CGALRS_dyadic_ptr rop,mpz_srcptr z)55 inline void CGALRS_dyadic_set_z(CGALRS_dyadic_ptr rop,mpz_srcptr z){
56         size_t prec;
57         prec=mpz_sizeinbase(z,2)-(mpz_tstbit(z,0)?0:mpz_scan1(z,0));
58         CGALRS_dyadic_set_prec(rop,prec);
59         mpfr_set_z(rop,z,GMP_RNDN);
60         CGAL_assertion(!mpfr_cmp_z(rop,z));
61 }
62 
CGALRS_dyadic_set_si(CGALRS_dyadic_ptr rop,long s)63 inline void CGALRS_dyadic_set_si(CGALRS_dyadic_ptr rop,long s){
64         CGALRS_dyadic_set_prec(rop,sizeof(long));
65         mpfr_set_si(rop,s,GMP_RNDN);
66         CGAL_assertion(!mpfr_cmp_si(rop,s));
67 }
68 
CGALRS_dyadic_set_ui(CGALRS_dyadic_ptr rop,unsigned long u)69 inline void CGALRS_dyadic_set_ui(CGALRS_dyadic_ptr rop,unsigned long u){
70         CGALRS_dyadic_set_prec(rop,sizeof(unsigned long));
71         mpfr_set_ui(rop,u,GMP_RNDN);
72         CGAL_assertion(!mpfr_cmp_ui(rop,u));
73 }
74 
CGALRS_dyadic_set_fr(CGALRS_dyadic_ptr rop,mpfr_srcptr op)75 inline void CGALRS_dyadic_set_fr(CGALRS_dyadic_ptr rop,mpfr_srcptr op){
76         if(rop!=op){
77                 CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
78                 mpfr_set(rop,op,GMP_RNDN);
79                 CGAL_assertion(mpfr_equal_p(rop,op)!=0);
80         }
81 }
82 
83 #define CGALRS_dyadic_init_set(R,D) \
84  ( CGALRS_dyadic_init(R), CGALRS_dyadic_set((R), (D)) )
85 #define CGALRS_dyadic_init_set_z(R,Z) \
86  ( CGALRS_dyadic_init(R), CGALRS_dyadic_set_z((R), (Z)) )
87 #define CGALRS_dyadic_init_set_si(R,I) \
88  ( CGALRS_dyadic_init(R), CGALRS_dyadic_set_si((R), (I)) )
89 #define CGALRS_dyadic_init_set_ui(R,I) \
90  ( CGALRS_dyadic_init(R), CGALRS_dyadic_set_ui((R), (I)) )
91 #define CGALRS_dyadic_init_set_fr(R,F) \
92  ( CGALRS_dyadic_init(R), CGALRS_dyadic_set_fr((R), (F)) )
93 
94 #define CGALRS_dyadic_get_fr(M,D)      mpfr_set(M,D,GMP_RNDN)
95 #define CGALRS_dyadic_get_d(D,RM)      mpfr_get_d(D,RM)
CGALRS_dyadic_get_exactfr(mpfr_ptr rop,CGALRS_dyadic_srcptr op)96 inline void CGALRS_dyadic_get_exactfr(mpfr_ptr rop,CGALRS_dyadic_srcptr op){
97         if(rop!=op){
98                 CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
99                 mpfr_set(rop,op,GMP_RNDN);
100                 CGAL_assertion(mpfr_equal_p(rop,op)!=0);
101         }
102 }
103 
104 #define CGALRS_dyadic_canonicalize(D)  ()
105 
106 // comparison functions
107 #define CGALRS_dyadic_sgn(D)    mpfr_sgn(D)
108 #define CGALRS_dyadic_zero(D)   mpfr_zero_p(D)
109 #define CGALRS_dyadic_cmp(D,E)  mpfr_cmp(D,E)
110 
111 // arithmetic functions
112 #define CGALRS_dyadic_add(R,D,E)        CGALRS_dyadic_ll_add(R,D,E,0)
113 #define CGALRS_dyadic_sub(R,D,E)        CGALRS_dyadic_ll_sub(R,D,E,0)
114 #define CGALRS_dyadic_mul(R,D,E)        CGALRS_dyadic_ll_mul(R,D,E,0)
115 
CGALRS_dyadic_neg(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op)116 inline void CGALRS_dyadic_neg(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
117         if(rop!=op)
118                 CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
119         mpfr_neg(rop,op,GMP_RNDN);
120         CGAL_assertion(
121                 rop==op||
122                 (!mpfr_cmpabs(rop,op)&&
123                 ((CGALRS_dyadic_zero(op)&&CGALRS_dyadic_zero(rop))||
124                  (CGALRS_dyadic_sgn(op)!=CGALRS_dyadic_sgn(rop)))));
125 }
126 
127 // low-level addition:
128 // add op1 and op2 and reserve b bits for future lowlevel operations
CGALRS_dyadic_ll_add(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,CGALRS_dyadic_srcptr op2,mp_prec_t b)129 inline void CGALRS_dyadic_ll_add(CGALRS_dyadic_ptr rop,
130                                  CGALRS_dyadic_srcptr op1,
131                                  CGALRS_dyadic_srcptr op2,
132                                  mp_prec_t b){
133         mp_exp_t l,r,temp1,temp2;
134         mp_prec_t rop_prec;
135         if(mpfr_zero_p(op1)){
136                 if(rop!=op2)
137                         CGALRS_dyadic_set(rop,op2);
138                 return;
139         }
140         if(mpfr_zero_p(op2)){
141                 if(rop!=op1)
142                         CGALRS_dyadic_set(rop,op1);
143                 return;
144         }
145         l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
146                 mpfr_get_exp(op1):
147                 mpfr_get_exp(op2);
148         temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
149         temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
150         r=temp1>temp2?temp2:temp1;
151         CGAL_assertion(l>r);
152 
153         rop_prec=b+1+(mp_prec_t)(l-r);
154         CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
155                 rop_prec>=mpfr_get_prec(op2));
156         if(rop==op1||rop==op2)
157                 CGALRS_dyadic_prec_round(rop,rop_prec);
158         else
159                 CGALRS_dyadic_set_prec(rop,rop_prec);
160         CGAL_assertion_code(int round=)
161         mpfr_add(rop,op1,op2,GMP_RNDN);
162         CGAL_assertion(!round);
163 }
164 
CGALRS_dyadic_add_z(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,mpz_srcptr z)165 inline void CGALRS_dyadic_add_z(CGALRS_dyadic_ptr rop,
166                                 CGALRS_dyadic_srcptr op1,
167                                 mpz_srcptr z){
168         mp_exp_t l,r;
169         mp_prec_t rop_prec;
170         if(mpfr_zero_p(op1)){
171                 CGALRS_dyadic_set_z(rop,z);
172                 return;
173         }
174         if(!mpz_sgn(z)){
175                 if(rop!=op1)
176                         CGALRS_dyadic_set(rop,op1);
177                 return;
178         }
179         l=mpfr_get_exp(op1)>(mp_exp_t)mpz_sizeinbase(z,2)?
180                 mpfr_get_exp(op1):
181                 mpz_sizeinbase(z,2);
182         r=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1)<0?
183                 mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1):
184                 0;
185         CGAL_assertion(l>r);
186 
187         rop_prec=1+(mp_prec_t)(l-r);
188         CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
189                 rop_prec>=(mp_prec_t)mpz_sizeinbase(z,2));
190         if(rop==op1)
191                 CGALRS_dyadic_prec_round(rop,rop_prec);
192         else
193                 CGALRS_dyadic_set_prec(rop,rop_prec);
194         CGAL_assertion_code(int round=)
195         mpfr_add_z(rop,op1,z,GMP_RNDN);
196         CGAL_assertion(!round);
197 }
198 
199 // low-level subtraction:
200 // subtract op2 to op1 and reserve b bits for future lowlevel operations
CGALRS_dyadic_ll_sub(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,CGALRS_dyadic_srcptr op2,mp_prec_t b)201 inline void CGALRS_dyadic_ll_sub(CGALRS_dyadic_ptr rop,
202                                  CGALRS_dyadic_srcptr op1,
203                                  CGALRS_dyadic_srcptr op2,
204                                  mp_prec_t b){
205         mp_exp_t l,r,temp1,temp2;
206         mp_prec_t rop_prec;
207         if(mpfr_zero_p(op1)){
208                 CGALRS_dyadic_neg(rop,op2);
209                 return;
210         }
211         if(mpfr_zero_p(op2)){
212                 if(rop!=op1)
213                         CGALRS_dyadic_set(rop,op1);
214                 return;
215         }
216         l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
217                 mpfr_get_exp(op1):
218                 mpfr_get_exp(op2);
219         temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
220         temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
221         r=temp1>temp2?temp2:temp1;
222         CGAL_assertion(l>r);
223 
224         rop_prec=b+1+(mp_prec_t)(l-r);
225         CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
226                 rop_prec>=mpfr_get_prec(op2));
227         if(rop==op1||rop==op2)
228                 CGALRS_dyadic_prec_round(rop,rop_prec);
229         else
230                 CGALRS_dyadic_set_prec(rop,rop_prec);
231         CGAL_assertion_code(int round=)
232         mpfr_sub(rop,op1,op2,GMP_RNDN);
233         CGAL_assertion(!round);
234 }
235 
236 // low-level multiplication:
237 // multiply op1 and op2 and reserve b bits for future lowlevel operations
CGALRS_dyadic_ll_mul(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,CGALRS_dyadic_srcptr op2,mp_prec_t b)238 inline void CGALRS_dyadic_ll_mul(CGALRS_dyadic_ptr rop,
239                                  CGALRS_dyadic_srcptr op1,
240                                  CGALRS_dyadic_srcptr op2,
241                                  mp_prec_t b){
242         if(rop==op1||rop==op2)
243                 CGALRS_dyadic_prec_round(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
244         else
245                 CGALRS_dyadic_set_prec(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
246         CGAL_assertion_code(int round=)
247         mpfr_mul(rop,op1,op2,GMP_RNDN);
248         CGAL_assertion(!round);
249 }
250 
CGALRS_dyadic_mul_z(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,mpz_srcptr z)251 inline void CGALRS_dyadic_mul_z(CGALRS_dyadic_ptr rop,
252                                 CGALRS_dyadic_srcptr op1,
253                                 mpz_srcptr z){
254         if(rop==op1)
255                 CGALRS_dyadic_prec_round(
256                         rop,
257                         mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
258         else
259                 CGALRS_dyadic_set_prec(
260                         rop,
261                         mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
262         CGAL_assertion_code(int round=)
263         mpfr_mul_z(rop,op1,z,GMP_RNDN);
264         CGAL_assertion(!round);
265 }
266 
CGALRS_dyadic_mul_si(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,long s)267 inline void CGALRS_dyadic_mul_si(CGALRS_dyadic_ptr rop,
268                                  CGALRS_dyadic_srcptr op1,
269                                  long s){
270         if(rop==op1)
271                 CGALRS_dyadic_prec_round(rop,mpfr_get_prec(op1)+sizeof(long));
272         else
273                 CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op1)+sizeof(long));
274         CGAL_assertion_code(int round=)
275         mpfr_mul_si(rop,op1,s,GMP_RNDN);
276         CGAL_assertion(!round);
277 }
278 
CGALRS_dyadic_mul_ui(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,unsigned long u)279 inline void CGALRS_dyadic_mul_ui(CGALRS_dyadic_ptr rop,
280                                  CGALRS_dyadic_srcptr op1,
281                                  unsigned long u){
282         if(rop==op1)
283                 CGALRS_dyadic_prec_round(
284                         rop,
285                         mpfr_get_prec(op1)+sizeof(unsigned long));
286         else
287                 CGALRS_dyadic_set_prec(
288                         rop,
289                         mpfr_get_prec(op1)+sizeof(unsigned long));
290         CGAL_assertion_code(int round=)
291         mpfr_mul_ui(rop,op1,u,GMP_RNDN);
292         CGAL_assertion(!round);
293 }
294 
CGALRS_dyadic_pow_ui(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,unsigned long u)295 inline void CGALRS_dyadic_pow_ui(CGALRS_dyadic_ptr rop,
296                                  CGALRS_dyadic_srcptr op1,
297                                  unsigned long u){
298         if(!u){
299                 CGAL_assertion_msg(!mpfr_zero_p(op1),"0^0");
300                 CGALRS_dyadic_set_ui(rop,1);
301                 return;
302         }
303         if(u==1){
304                 if(rop!=op1)
305                         CGALRS_dyadic_set(rop,op1);
306                 return;
307         }
308         if(mpfr_zero_p(op1)){
309                 CGAL_assertion_msg(u!=0,"0^0");
310                 CGALRS_dyadic_set_ui(rop,0);
311                 return;
312         }
313         if(!mpfr_cmp_ui(op1,1)){
314                 if(rop!=op1)
315                         CGALRS_dyadic_set(rop,op1);
316                 return;
317         }
318         if(rop==op1)
319                 CGALRS_dyadic_prec_round(rop,u*mpfr_get_prec(op1));
320         else
321                 CGALRS_dyadic_set_prec(rop,u*mpfr_get_prec(op1));
322         CGAL_assertion_code(int round=)
323         mpfr_pow_ui(rop,op1,u,GMP_RNDN);
324         CGAL_assertion(!round);
325 }
326 
CGALRS_dyadic_addmul(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,CGALRS_dyadic_srcptr op2)327 inline void CGALRS_dyadic_addmul(CGALRS_dyadic_ptr rop,
328                                  CGALRS_dyadic_srcptr op1,
329                                  CGALRS_dyadic_srcptr op2){
330         CGALRS_dyadic_t temp;
331         CGALRS_dyadic_init(temp);
332         CGALRS_dyadic_mul(temp,op1,op2);
333         CGALRS_dyadic_add(rop,rop,temp);
334         CGALRS_dyadic_clear(temp);
335 }
336 
CGALRS_dyadic_addmul_si(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,long op2)337 inline void CGALRS_dyadic_addmul_si(CGALRS_dyadic_ptr rop,
338                                     CGALRS_dyadic_srcptr op1,
339                                     long op2){
340         CGALRS_dyadic_t temp;
341         CGALRS_dyadic_init(temp);
342         CGALRS_dyadic_mul_si(temp,op1,op2);
343         CGALRS_dyadic_add(rop,rop,temp);
344         CGALRS_dyadic_clear(temp);
345 }
346 
CGALRS_dyadic_addmul_ui(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,unsigned long u)347 inline void CGALRS_dyadic_addmul_ui(CGALRS_dyadic_ptr rop,
348                                     CGALRS_dyadic_srcptr op1,
349                                     unsigned long u){
350         //CGALRS_dyadic_t temp;
351         //CGALRS_dyadic_init(temp);
352         //CGALRS_dyadic_mul_ui(temp,op1,u);
353         //CGALRS_dyadic_add(rop,rop,temp);
354         //CGALRS_dyadic_clear(temp);
355         CGALRS_dyadic_t temp;
356         mp_exp_t l,r,temp1,temp2;
357         mp_prec_t rop_prec;
358         if(u==0||mpfr_zero_p(op1))
359                 return;
360         if(u==1){
361                 CGALRS_dyadic_add(rop,rop,op1);
362                 return;
363         }
364         // TODO: if(op1==1)
365         // calculate temp=op1*u
366         mpfr_init2(temp,mpfr_get_prec(op1)+sizeof(unsigned int));
367         CGAL_assertion_code(int round1=)
368         mpfr_mul_ui(temp,op1,u,GMP_RNDN);
369         CGAL_assertion(!round1);
370         // calculate the precision needed for rop
371         l=mpfr_get_exp(op1)>0?mpfr_get_exp(op1):0;
372         temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
373         temp2=sizeof(unsigned long);
374         r=temp1>temp2?temp2:temp1;
375         CGAL_assertion(l>r);
376         rop_prec=sizeof(unsigned long)+1+(mp_prec_t)(l-r);
377         CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
378                 rop_prec>=mpfr_get_prec(rop));
379         // set precision and add
380         CGALRS_dyadic_prec_round(rop,rop_prec);
381         CGAL_assertion_code(int round2=)
382         mpfr_add(rop,rop,temp,GMP_RNDN);
383         CGAL_assertion(!round2);
384 }
385 
CGALRS_dyadic_mul_2exp(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,unsigned long ui)386 inline void CGALRS_dyadic_mul_2exp(CGALRS_dyadic_ptr rop,
387                                    CGALRS_dyadic_srcptr op1,
388                                    unsigned long ui){
389         // mpfr_mul_2ui does change the mantissa!
390         if(rop==op1)
391                 CGALRS_dyadic_prec_round(
392                         rop,
393                         sizeof(unsigned long)+mpfr_get_prec(op1));
394         else
395                 CGALRS_dyadic_set_prec(
396                         rop,
397                         sizeof(unsigned long)+mpfr_get_prec(op1));
398         CGAL_assertion_code(int round=)
399         mpfr_mul_2ui(rop,op1,ui,GMP_RNDN);
400         CGAL_assertion(!round);
401 }
402 
CGALRS_dyadic_div_2exp(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op1,unsigned long ui)403 inline void CGALRS_dyadic_div_2exp(CGALRS_dyadic_ptr rop,
404                                    CGALRS_dyadic_srcptr op1,
405                                    unsigned long ui){
406         // mpfr_div_2ui does not change the mantissa... am I sure?
407         CGAL_assertion_code(int round=)
408         mpfr_div_2ui(rop,op1,ui,GMP_RNDN);
409         CGAL_assertion(!round);
410 }
411 
412 // miscellaneous functions
413 #define CGALRS_dyadic_midpoint(R,D,E) \
414  ( CGALRS_dyadic_ll_add(R,D,E,1) , mpfr_div_2ui(R,R,1,GMP_RNDN) )
415 #define CGALRS_dyadic_swap(D,E)         mpfr_swap(D,E)
416 
417 // I/O functions
418 #define CGALRS_dyadic_out_str(F,D)      mpfr_out_str(F,10,0,D,GMP_RNDN)
419 #ifdef __cplusplus
420 inline std::ostream& operator<<(std::ostream &s,CGALRS_dyadic_srcptr op){
421         mp_exp_t exponent;
422         mpz_t mantissa;
423         mpz_init(mantissa);
424         exponent=mpfr_get_z_exp(mantissa,op);
425         s<<"["<<mantissa<<","<<exponent<<"]";
426         return s;
427 }
428 #endif
429 
430 #endif  // CGAL_RS_DYADIC_H
431