1 /* mini-mpq, a minimalistic implementation of a GNU GMP subset.
2
3 Contributed to the GNU project by Marco Bodrato
4
5 Acknowledgment: special thanks to Bradley Lucier for his comments
6 to the preliminary version of this code.
7
8 Copyright 2018-2022 Free Software Foundation, Inc.
9
10 This file is part of the GNU MP Library.
11
12 The GNU MP Library is free software; you can redistribute it and/or modify
13 it under the terms of either:
14
15 * the GNU Lesser General Public License as published by the Free
16 Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
18
19 or
20
21 * the GNU General Public License as published by the Free Software
22 Foundation; either version 2 of the License, or (at your option) any
23 later version.
24
25 or both in parallel, as here.
26
27 The GNU MP Library is distributed in the hope that it will be useful, but
28 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
30 for more details.
31
32 You should have received copies of the GNU General Public License and the
33 GNU Lesser General Public License along with the GNU MP Library. If not,
34 see https://www.gnu.org/licenses/. */
35
36 #include <assert.h>
37 #include <limits.h>
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42 #include "mini-mpq.h"
43
44 #ifndef GMP_LIMB_HIGHBIT
45 /* Define macros and static functions already defined by mini-gmp.c */
46 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
47 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
48 #define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0)
49 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
50 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
51
52 static mpz_srcptr
53 mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
54 {
55 x->_mp_alloc = 0;
56 x->_mp_d = (mp_ptr) xp;
57 x->_mp_size = xs;
58 return x;
59 }
60
61 static void
62 gmp_die (const char *msg)
63 {
64 fprintf (stderr, "%s\n", msg);
65 abort();
66 }
67 #endif
68
69
70 /* MPQ helper functions */
71 static mpq_srcptr
72 mpq_roinit_normal_nn (mpq_t x, mp_srcptr np, mp_size_t ns,
73 mp_srcptr dp, mp_size_t ds)
74 {
75 mpz_roinit_normal_n (mpq_numref(x), np, ns);
76 mpz_roinit_normal_n (mpq_denref(x), dp, ds);
77 return x;
78 }
79
80 static mpq_srcptr
81 mpq_roinit_zz (mpq_t x, mpz_srcptr n, mpz_srcptr d)
82 {
83 return mpq_roinit_normal_nn (x, n->_mp_d, n->_mp_size,
84 d->_mp_d, d->_mp_size);
85 }
86
87 static void
88 mpq_nan_init (mpq_t x)
89 {
90 mpz_init (mpq_numref (x));
91 mpz_init (mpq_denref (x));
92 }
93
94 void
95 mpq_init (mpq_t x)
96 {
97 mpz_init (mpq_numref (x));
98 mpz_init_set_ui (mpq_denref (x), 1);
99 }
100
101 void
102 mpq_clear (mpq_t x)
103 {
104 mpz_clear (mpq_numref (x));
105 mpz_clear (mpq_denref (x));
106 }
107
108 static void
109 mpq_canonical_sign (mpq_t r)
110 {
111 mp_size_t ds = mpq_denref (r)->_mp_size;
112 if (ds <= 0)
113 {
114 if (ds == 0)
115 gmp_die("mpq: Fraction with zero denominator.");
116 mpz_neg (mpq_denref (r), mpq_denref (r));
117 mpz_neg (mpq_numref (r), mpq_numref (r));
118 }
119 }
120
121 static void
122 mpq_helper_canonicalize (mpq_t r, const mpz_t num, const mpz_t den)
123 {
124 if (num->_mp_size == 0)
125 mpq_set_ui (r, 0, 1);
126 else
127 {
128 mpz_t g;
129
130 mpz_init (g);
131 mpz_gcd (g, num, den);
132 mpz_tdiv_q (mpq_numref (r), num, g);
133 mpz_tdiv_q (mpq_denref (r), den, g);
134 mpz_clear (g);
135 mpq_canonical_sign (r);
136 }
137 }
138
139 void
140 mpq_canonicalize (mpq_t r)
141 {
142 mpq_helper_canonicalize (r, mpq_numref (r), mpq_denref (r));
143 }
144
145 void
146 mpq_swap (mpq_t a, mpq_t b)
147 {
148 mpz_swap (mpq_numref (a), mpq_numref (b));
149 mpz_swap (mpq_denref (a), mpq_denref (b));
150 }
151
152
153 /* MPQ assignment and conversions. */
154 void
155 mpz_set_q (mpz_t r, const mpq_t q)
156 {
157 mpz_tdiv_q (r, mpq_numref (q), mpq_denref (q));
158 }
159
160 void
161 mpq_set (mpq_t r, const mpq_t q)
162 {
163 mpz_set (mpq_numref (r), mpq_numref (q));
164 mpz_set (mpq_denref (r), mpq_denref (q));
165 }
166
167 void
168 mpq_set_ui (mpq_t r, unsigned long n, unsigned long d)
169 {
170 mpz_set_ui (mpq_numref (r), n);
171 mpz_set_ui (mpq_denref (r), d);
172 }
173
174 void
175 mpq_set_si (mpq_t r, signed long n, unsigned long d)
176 {
177 mpz_set_si (mpq_numref (r), n);
178 mpz_set_ui (mpq_denref (r), d);
179 }
180
181 void
182 mpq_set_z (mpq_t r, const mpz_t n)
183 {
184 mpz_set_ui (mpq_denref (r), 1);
185 mpz_set (mpq_numref (r), n);
186 }
187
188 void
189 mpq_set_num (mpq_t r, const mpz_t z)
190 {
191 mpz_set (mpq_numref (r), z);
192 }
193
194 void
195 mpq_set_den (mpq_t r, const mpz_t z)
196 {
197 mpz_set (mpq_denref (r), z);
198 }
199
200 void
201 mpq_get_num (mpz_t r, const mpq_t q)
202 {
203 mpz_set (r, mpq_numref (q));
204 }
205
206 void
207 mpq_get_den (mpz_t r, const mpq_t q)
208 {
209 mpz_set (r, mpq_denref (q));
210 }
211
212
213 /* MPQ comparisons and the like. */
214 int
215 mpq_cmp (const mpq_t a, const mpq_t b)
216 {
217 mpz_t t1, t2;
218 int res;
219
220 mpz_init (t1);
221 mpz_init (t2);
222 mpz_mul (t1, mpq_numref (a), mpq_denref (b));
223 mpz_mul (t2, mpq_numref (b), mpq_denref (a));
224 res = mpz_cmp (t1, t2);
225 mpz_clear (t1);
226 mpz_clear (t2);
227
228 return res;
229 }
230
231 int
232 mpq_cmp_z (const mpq_t a, const mpz_t b)
233 {
234 mpz_t t;
235 int res;
236
237 mpz_init (t);
238 mpz_mul (t, b, mpq_denref (a));
239 res = mpz_cmp (mpq_numref (a), t);
240 mpz_clear (t);
241
242 return res;
243 }
244
245 int
246 mpq_equal (const mpq_t a, const mpq_t b)
247 {
248 return (mpz_cmp (mpq_numref (a), mpq_numref (b)) == 0) &&
249 (mpz_cmp (mpq_denref (a), mpq_denref (b)) == 0);
250 }
251
252 int
253 mpq_cmp_ui (const mpq_t q, unsigned long n, unsigned long d)
254 {
255 mpq_t t;
256 assert (d != 0);
257 if (ULONG_MAX <= GMP_LIMB_MAX) {
258 mp_limb_t nl = n, dl = d;
259 return mpq_cmp (q, mpq_roinit_normal_nn (t, &nl, n != 0, &dl, 1));
260 } else {
261 int ret;
262
263 mpq_nan_init (t);
264 mpq_set_ui (t, n, d);
265 ret = mpq_cmp (q, t);
266 mpq_clear (t);
267
268 return ret;
269 }
270 }
271
272 int
273 mpq_cmp_si (const mpq_t q, signed long n, unsigned long d)
274 {
275 assert (d != 0);
276
277 if (n >= 0)
278 return mpq_cmp_ui (q, n, d);
279 else
280 {
281 mpq_t t;
282
283 if (ULONG_MAX <= GMP_LIMB_MAX)
284 {
285 mp_limb_t nl = GMP_NEG_CAST (unsigned long, n), dl = d;
286 return mpq_cmp (q, mpq_roinit_normal_nn (t, &nl, -1, &dl, 1));
287 }
288 else
289 {
290 unsigned long l_n = GMP_NEG_CAST (unsigned long, n);
291
292 mpq_roinit_normal_nn (t, mpq_numref (q)->_mp_d, - mpq_numref (q)->_mp_size,
293 mpq_denref (q)->_mp_d, mpq_denref (q)->_mp_size);
294 return - mpq_cmp_ui (t, l_n, d);
295 }
296 }
297 }
298
299 int
300 mpq_sgn (const mpq_t a)
301 {
302 return mpz_sgn (mpq_numref (a));
303 }
304
305
306 /* MPQ arithmetic. */
307 void
308 mpq_abs (mpq_t r, const mpq_t q)
309 {
310 mpz_abs (mpq_numref (r), mpq_numref (q));
311 mpz_set (mpq_denref (r), mpq_denref (q));
312 }
313
314 void
315 mpq_neg (mpq_t r, const mpq_t q)
316 {
317 mpz_neg (mpq_numref (r), mpq_numref (q));
318 mpz_set (mpq_denref (r), mpq_denref (q));
319 }
320
321 void
322 mpq_add (mpq_t r, const mpq_t a, const mpq_t b)
323 {
324 mpz_t t;
325
326 mpz_init (t);
327 mpz_gcd (t, mpq_denref (a), mpq_denref (b));
328 if (mpz_cmp_ui (t, 1) == 0)
329 {
330 mpz_mul (t, mpq_numref (a), mpq_denref (b));
331 mpz_addmul (t, mpq_numref (b), mpq_denref (a));
332 mpz_mul (mpq_denref (r), mpq_denref (a), mpq_denref (b));
333 mpz_swap (mpq_numref (r), t);
334 }
335 else
336 {
337 mpz_t x, y;
338 mpz_init (x);
339 mpz_init (y);
340
341 mpz_tdiv_q (x, mpq_denref (b), t);
342 mpz_tdiv_q (y, mpq_denref (a), t);
343 mpz_mul (x, mpq_numref (a), x);
344 mpz_addmul (x, mpq_numref (b), y);
345
346 mpz_gcd (t, x, t);
347 mpz_tdiv_q (mpq_numref (r), x, t);
348 mpz_tdiv_q (x, mpq_denref (b), t);
349 mpz_mul (mpq_denref (r), x, y);
350
351 mpz_clear (x);
352 mpz_clear (y);
353 }
354 mpz_clear (t);
355 }
356
357 void
358 mpq_sub (mpq_t r, const mpq_t a, const mpq_t b)
359 {
360 mpq_t t;
361
362 mpq_roinit_normal_nn (t, mpq_numref (b)->_mp_d, - mpq_numref (b)->_mp_size,
363 mpq_denref (b)->_mp_d, mpq_denref (b)->_mp_size);
364 mpq_add (r, a, t);
365 }
366
367 void
368 mpq_div (mpq_t r, const mpq_t a, const mpq_t b)
369 {
370 mpq_t t;
371 mpq_mul (r, a, mpq_roinit_zz (t, mpq_denref (b), mpq_numref (b)));
372 }
373
374 void
375 mpq_mul (mpq_t r, const mpq_t a, const mpq_t b)
376 {
377 mpq_t t;
378 mpq_nan_init (t);
379
380 if (a != b) {
381 mpq_helper_canonicalize (t, mpq_numref (a), mpq_denref (b));
382 mpq_helper_canonicalize (r, mpq_numref (b), mpq_denref (a));
383
384 a = r;
385 b = t;
386 }
387
388 mpz_mul (mpq_numref (r), mpq_numref (a), mpq_numref (b));
389 mpz_mul (mpq_denref (r), mpq_denref (a), mpq_denref (b));
390 mpq_clear (t);
391 }
392
393 static void
394 mpq_helper_2exp (mpz_t rn, mpz_t rd, const mpz_t qn, const mpz_t qd, mp_bitcnt_t e)
395 {
396 mp_bitcnt_t z = mpz_scan1 (qd, 0);
397 z = GMP_MIN (z, e);
398 mpz_mul_2exp (rn, qn, e - z);
399 mpz_tdiv_q_2exp (rd, qd, z);
400 }
401
402 void
403 mpq_div_2exp (mpq_t r, const mpq_t q, mp_bitcnt_t e)
404 {
405 mpq_helper_2exp (mpq_denref (r), mpq_numref (r), mpq_denref (q), mpq_numref (q), e);
406 }
407
408 void
409 mpq_mul_2exp (mpq_t r, const mpq_t q, mp_bitcnt_t e)
410 {
411 mpq_helper_2exp (mpq_numref (r), mpq_denref (r), mpq_numref (q), mpq_denref (q), e);
412 }
413
414 void
415 mpq_inv (mpq_t r, const mpq_t q)
416 {
417 mpq_set (r, q);
418 mpz_swap (mpq_denref (r), mpq_numref (r));
419 mpq_canonical_sign (r);
420 }
421
422
423 /* MPQ to/from double. */
424 void
425 mpq_set_d (mpq_t r, double x)
426 {
427 mpz_set_ui (mpq_denref (r), 1);
428
429 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
430 zero or infinity. */
431 if (x == x * 0.5 || x != x)
432 mpq_numref (r)->_mp_size = 0;
433 else
434 {
435 double B;
436 mp_bitcnt_t e;
437
438 B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
439 for (e = 0; x != x + 0.5; e += GMP_LIMB_BITS)
440 x *= B;
441
442 mpz_set_d (mpq_numref (r), x);
443 mpq_div_2exp (r, r, e);
444 }
445 }
446
447 double
448 mpq_get_d (const mpq_t u)
449 {
450 mp_bitcnt_t ne, de, ee;
451 mpz_t z;
452 double B, ret;
453
454 ne = mpz_sizeinbase (mpq_numref (u), 2);
455 de = mpz_sizeinbase (mpq_denref (u), 2);
456
457 ee = CHAR_BIT * sizeof (double);
458 if (de == 1 || ne > de + ee)
459 ee = 0;
460 else
461 ee = (ee + de - ne) / GMP_LIMB_BITS + 1;
462
463 mpz_init (z);
464 mpz_mul_2exp (z, mpq_numref (u), ee * GMP_LIMB_BITS);
465 mpz_tdiv_q (z, z, mpq_denref (u));
466 ret = mpz_get_d (z);
467 mpz_clear (z);
468
469 B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
470 for (B = 1 / B; ee != 0; --ee)
471 ret *= B;
472
473 return ret;
474 }
475
476
477 /* MPQ and strings/streams. */
478 char *
479 mpq_get_str (char *sp, int base, const mpq_t q)
480 {
481 char *res;
482 char *rden;
483 size_t len;
484
485 res = mpz_get_str (sp, base, mpq_numref (q));
486 if (res == NULL || mpz_cmp_ui (mpq_denref (q), 1) == 0)
487 return res;
488
489 len = strlen (res) + 1;
490 rden = sp ? sp + len : NULL;
491 rden = mpz_get_str (rden, base, mpq_denref (q));
492 assert (rden != NULL);
493
494 if (sp == NULL) {
495 void * (*gmp_reallocate_func) (void *, size_t, size_t);
496 void (*gmp_free_func) (void *, size_t);
497 size_t lden;
498
499 mp_get_memory_functions (NULL, &gmp_reallocate_func, &gmp_free_func);
500 lden = strlen (rden) + 1;
501 res = (char *) gmp_reallocate_func (res, len, (lden + len) * sizeof (char));
502 memcpy (res + len, rden, lden);
503 gmp_free_func (rden, lden);
504 }
505
506 res [len - 1] = '/';
507 return res;
508 }
509
510 size_t
511 mpq_out_str (FILE *stream, int base, const mpq_t x)
512 {
513 char * str;
514 size_t len, n;
515 void (*gmp_free_func) (void *, size_t);
516
517 str = mpq_get_str (NULL, base, x);
518 if (!str)
519 return 0;
520 len = strlen (str);
521 n = fwrite (str, 1, len, stream);
522 mp_get_memory_functions (NULL, NULL, &gmp_free_func);
523 gmp_free_func (str, len + 1);
524 return n;
525 }
526
527 int
528 mpq_set_str (mpq_t r, const char *sp, int base)
529 {
530 const char *slash;
531
532 slash = strchr (sp, '/');
533 if (slash == NULL) {
534 mpz_set_ui (mpq_denref(r), 1);
535 return mpz_set_str (mpq_numref(r), sp, base);
536 } else {
537 char *num;
538 size_t numlen;
539 int ret;
540 void * (*gmp_allocate_func) (size_t);
541 void (*gmp_free_func) (void *, size_t);
542
543 mp_get_memory_functions (&gmp_allocate_func, NULL, &gmp_free_func);
544 numlen = slash - sp;
545 num = (char *) gmp_allocate_func (numlen + 1);
546 memcpy (num, sp, numlen);
547 num[numlen] = '\0';
548 ret = mpz_set_str (mpq_numref(r), num, base);
549 gmp_free_func (num, numlen + 1);
550
551 if (ret != 0)
552 return ret;
553
554 return mpz_set_str (mpq_denref(r), slash + 1, base);
555 }
556 }