1 /*
2 * This file is part of MPSolve 3.2.1
3 *
4 * Copyright (C) 2001-2020, Dipartimento di Matematica "L. Tonelli", Pisa.
5 * License: http://www.gnu.org/licenses/gpl.html GPL version 3 or higher
6 *
7 * Authors:
8 * Dario Andrea Bini <bini@dm.unipi.it>
9 * Giuseppe Fiorentino <fiorent@dm.unipi.it>
10 * Leonardo Robol <leonardo.robol@unipi.it>
11 */
12
13
14 #include <stdlib.h>
15 #include <math.h>
16 #include <mps/mps.h>
17
18 /**********************************************
19 * MPZ_T *
20 **********************************************/
21
22 #ifndef mpz_swap
23 void
mpz_swap(mpz_t z1,mpz_t z2)24 mpz_swap (mpz_t z1, mpz_t z2)
25 /* z1 <-> z2 */
26 {
27 mpz_t t;
28
29 mpz_Move (t, z1);
30 mpz_Move (z1, z2);
31 mpz_Move (z2, t);
32 }
33 #endif
34
35 #ifndef mpz_tstbit
36 int
mpz_tstbit(mpz_t z,unsigned long int pos)37 mpz_tstbit (mpz_t z, unsigned long int pos)
38 {
39 size_t size;
40
41 size = mpz_sizeinbase (z, 2);
42
43 if (mpz_sgn (z) == 0 || size <= pos)
44 return 0;
45
46 if (mpz_scan1 (z, pos) == pos)
47 return 1;
48 else
49 return 0;
50 }
51 #endif
52
53 /* vector support functions */
54
55 void
mpz_vinit(mpz_t v[],unsigned long int size)56 mpz_vinit (mpz_t v[], unsigned long int size)
57 {
58 unsigned long int i;
59
60 for (i = 0; i < size; i++)
61 mpz_init (v[i]);
62 }
63
64 void
mpz_vclear(mpz_t v[],unsigned long int size)65 mpz_vclear (mpz_t v[], unsigned long int size)
66 {
67 unsigned long int i;
68
69 for (i = 0; i < size; i++)
70 mpz_clear (v[i]);
71 }
72
73 /**********************************************
74 * MPQ_T *
75 **********************************************/
76
77 #ifndef mpq_swap
78 void
mpq_swap(mpq_t q1,mpq_t q2)79 mpq_swap (mpq_t q1, mpq_t q2)
80 /* q1 <-> q2 */
81 {
82 mpq_t t;
83
84 mpq_Move (t, q1);
85 mpq_Move (q1, q2);
86 mpq_Move (q2, t);
87 }
88 #endif
89
90 #ifndef mpq_out_str
91 void
mpq_out_str(FILE * stream,int base,mpq_t q)92 mpq_out_str (FILE * stream, int base, mpq_t q)
93 {
94 mpz_out_str (stdout, base, mpq_numref (q));
95 putc ('/', stream);
96 mpz_out_str (stdout, base, mpq_denref (q));
97 }
98 #endif
99
100 /* vector support functions */
101
102 void
mpq_vinit(mpq_t v[],unsigned long int size)103 mpq_vinit (mpq_t v[], unsigned long int size)
104 {
105 unsigned long int i;
106
107 for (i = 0; i < size; i++)
108 mpq_init (v[i]);
109 }
110
111 void
mpq_vclear(mpq_t v[],unsigned long int size)112 mpq_vclear (mpq_t v[], unsigned long int size)
113 {
114 unsigned long int i;
115
116 for (i = 0; i < size; i++)
117 mpq_clear (v[i]);
118 }
119
120 /**********************************************
121 * MPF_T *
122 **********************************************/
123 #ifndef mpf_swap
124 void
mpf_swap(mpf_t f1,mpf_t f2)125 mpf_swap (mpf_t f1, mpf_t f2)
126 /* f1 <-> f2 */
127 {
128 mpf_t t;
129
130 mpf_Move (t, f1);
131 mpf_Move (f1, f2);
132 mpf_Move (f2, t);
133 }
134 #endif
135
136 void
mpf_set_2dl(mpf_t f,double d,long int l)137 mpf_set_2dl (mpf_t f, double d, long int l)
138 {
139 mpf_set_d (f, d);
140 if (l >= 0)
141 mpf_mul_2exp (f, f, l);
142 else
143 mpf_div_2exp (f, f, -l);
144 }
145
146 void
mpf_get_2dl(double * d,long int * l,mpf_t f)147 mpf_get_2dl (double *d, long int *l, mpf_t f)
148 {
149 mp_exp_t e;
150 double t;
151 int i;
152
153 /* pick mantissa and exponent from f */
154 e = f->_mp_exp;
155 f->_mp_exp = 0;
156 t = mpf_get_d (f);
157 f->_mp_exp = e;
158
159 /* scale mantissa to (1/2, 1] */
160 *d = frexp (t, &i);
161 *l = e * mp_bits_per_limb + i;
162 }
163
164 long int
mpf_size_2(mpf_t f)165 mpf_size_2 (mpf_t f)
166 {
167 double d;
168 long int l;
169
170 /* pick mantissa and exponent from f */
171 mpf_get_2dl (&d, &l, f);
172
173 return l;
174 }
175
176 /* missing functions */
177
178 void
mpf_add_si(mpf_t r,mpf_t f,long int i)179 mpf_add_si (mpf_t r, mpf_t f, long int i)
180 {
181 if (i >= 0)
182 mpf_add_ui (r, f, i);
183 else
184 mpf_sub_ui (r, f, -i);
185 }
186
187 void
mpf_sub_si(mpf_t r,mpf_t f,long int i)188 mpf_sub_si (mpf_t r, mpf_t f, long int i)
189 {
190 if (i >= 0)
191 mpf_sub_ui (r, f, i);
192 else
193 mpf_add_ui (r, f, -i);
194 }
195
196 void
mpf_si_sub(mpf_t r,long int i,mpf_t f)197 mpf_si_sub (mpf_t r, long int i, mpf_t f)
198 {
199 if (i >= 0)
200 mpf_ui_sub (r, i, f);
201 else
202 {
203 mpf_add_ui (r, f, -i);
204 mpf_neg (r, r);
205 }
206 }
207
208 void
mpf_mul_si(mpf_t r,mpf_t f,long int i)209 mpf_mul_si (mpf_t r, mpf_t f, long int i)
210 {
211 if (i >= 0)
212 mpf_mul_ui (r, f, i);
213 else
214 {
215 mpf_mul_ui (r, f, -i);
216 mpf_neg (r, r);
217 }
218 }
219
220 void
mpf_div_si(mpf_t r,mpf_t f,long int i)221 mpf_div_si (mpf_t r, mpf_t f, long int i)
222 {
223 if (i >= 0)
224 mpf_div_ui (r, f, i);
225 else
226 {
227 mpf_div_ui (r, f, -i);
228 mpf_neg (r, r);
229 }
230 }
231
232 void
mpf_si_div(mpf_t r,long int i,mpf_t f)233 mpf_si_div (mpf_t r, long int i, mpf_t f)
234 {
235 if (i >= 0)
236 mpf_ui_div (r, i, f);
237 else
238 {
239 mpf_ui_div (r, -i, f);
240 mpf_neg (r, r);
241 }
242 }
243
244 #ifndef mpf_pow_ui
245 void
mpf_pow_ui(mpf_t r,mpf_t f,register unsigned long int i)246 mpf_pow_ui (mpf_t r, mpf_t f, register unsigned long int i)
247 /* r = f^i */
248 {
249 mpf_t t;
250
251 mpf_init2 (t, mpf_get_prec (r));
252
253 mpf_set (t, f);
254
255 if (i & 1)
256 mpf_set (r, t);
257 else
258 mpf_set_ui (r, 1);
259 i >>= 1; /* divide i by 2 */
260
261 while (i)
262 {
263 mpf_sqr_eq (t);
264 if (i & 1)
265 mpf_mul_eq (r, t);
266 i >>= 1; /* divide i by 2 */
267 }
268
269 mpf_clear (t);
270 }
271 #endif
272
273 void
mpf_pow_si(mpf_t r,mpf_t f,long int i)274 mpf_pow_si (mpf_t r, mpf_t f, long int i)
275 /* r = f^i, i integer */
276 {
277 if (i >= 0)
278 mpf_pow_ui (r, f, i);
279 else
280 {
281 mpf_pow_ui (r, f, -i);
282 mpf_inv_eq (r);
283 }
284 }
285
286 /* vector support functions */
287
288 void
mpf_vinit(mpf_t * v,unsigned long int size)289 mpf_vinit (mpf_t * v, unsigned long int size)
290 {
291 unsigned long int i;
292
293 for (i = 0; i < size; i++)
294 mpf_init (v[i]);
295 }
296
297 void
mpf_vinit2(mpf_t v[],unsigned long int size,unsigned long int prec)298 mpf_vinit2 (mpf_t v[], unsigned long int size, unsigned long int prec)
299 {
300 unsigned long int i;
301
302 for (i = 0; i < size; i++)
303 mpf_init2 (v[i], prec);
304 }
305
306 void
mpf_vclear(mpf_t v[],unsigned long int size)307 mpf_vclear (mpf_t v[], unsigned long int size)
308 {
309 unsigned long int i;
310
311 for (i = 0; i < size; i++)
312 mpf_clear (v[i]);
313 }
314