1/*
2 Name: imrat.c
3 Purpose: Arbitrary precision rational arithmetic routines.
4 Author: M. J. Fromberger
5
6 Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this software and associated documentation files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24 SOFTWARE.
25 */
26
27#include "imrat.h"
28#include <assert.h>
29#include <ctype.h>
30#include <stdlib.h>
31#include <string.h>
32
33#define MP_NUMER_SIGN(Q) (MP_NUMER_P(Q)->sign)
34#define MP_DENOM_SIGN(Q) (MP_DENOM_P(Q)->sign)
35
36#define TEMP(K) (temp + (K))
37#define SETUP(E, C) \
38 do { \
39 if ((res = (E)) != MP_OK) goto CLEANUP; \
40 ++(C); \
41 } while (0)
42
43/* Reduce the given rational, in place, to lowest terms and canonical form.
44 Zero is represented as 0/1, one as 1/1. Signs are adjusted so that the sign
45 of the numerator is definitive. */
46static mp_result s_rat_reduce(mp_rat r);
47
48/* Common code for addition and subtraction operations on rationals. */
49static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,
50 mp_result (*comb_f)(mp_int, mp_int, mp_int));
51
52mp_result mp_rat_init(mp_rat r) {
53 mp_result res;
54
55 if ((res = mp_int_init(z: MP_NUMER_P(Q: r))) != MP_OK) return res;
56 if ((res = mp_int_init(z: MP_DENOM_P(Q: r))) != MP_OK) {
57 mp_int_clear(z: MP_NUMER_P(Q: r));
58 return res;
59 }
60
61 return mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
62}
63
64mp_rat mp_rat_alloc(void) {
65 mp_rat out = malloc(size: sizeof(*out));
66
67 if (out != NULL) {
68 if (mp_rat_init(r: out) != MP_OK) {
69 free(ptr: out);
70 return NULL;
71 }
72 }
73
74 return out;
75}
76
77mp_result mp_rat_reduce(mp_rat r) { return s_rat_reduce(r); }
78
79mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec) {
80 mp_result res;
81
82 if ((res = mp_int_init_size(z: MP_NUMER_P(Q: r), prec: n_prec)) != MP_OK) {
83 return res;
84 }
85 if ((res = mp_int_init_size(z: MP_DENOM_P(Q: r), prec: d_prec)) != MP_OK) {
86 mp_int_clear(z: MP_NUMER_P(Q: r));
87 return res;
88 }
89
90 return mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
91}
92
93mp_result mp_rat_init_copy(mp_rat r, mp_rat old) {
94 mp_result res;
95
96 if ((res = mp_int_init_copy(z: MP_NUMER_P(Q: r), old: MP_NUMER_P(Q: old))) != MP_OK) {
97 return res;
98 }
99 if ((res = mp_int_init_copy(z: MP_DENOM_P(Q: r), old: MP_DENOM_P(Q: old))) != MP_OK)
100 mp_int_clear(z: MP_NUMER_P(Q: r));
101
102 return res;
103}
104
105mp_result mp_rat_set_value(mp_rat r, mp_small numer, mp_small denom) {
106 mp_result res;
107
108 if (denom == 0) return MP_UNDEF;
109
110 if ((res = mp_int_set_value(z: MP_NUMER_P(Q: r), value: numer)) != MP_OK) {
111 return res;
112 }
113 if ((res = mp_int_set_value(z: MP_DENOM_P(Q: r), value: denom)) != MP_OK) {
114 return res;
115 }
116
117 return s_rat_reduce(r);
118}
119
120mp_result mp_rat_set_uvalue(mp_rat r, mp_usmall numer, mp_usmall denom) {
121 mp_result res;
122
123 if (denom == 0) return MP_UNDEF;
124
125 if ((res = mp_int_set_uvalue(z: MP_NUMER_P(Q: r), uvalue: numer)) != MP_OK) {
126 return res;
127 }
128 if ((res = mp_int_set_uvalue(z: MP_DENOM_P(Q: r), uvalue: denom)) != MP_OK) {
129 return res;
130 }
131
132 return s_rat_reduce(r);
133}
134
135void mp_rat_clear(mp_rat r) {
136 mp_int_clear(z: MP_NUMER_P(Q: r));
137 mp_int_clear(z: MP_DENOM_P(Q: r));
138}
139
140void mp_rat_free(mp_rat r) {
141 assert(r != NULL);
142
143 if (r->num.digits != NULL) mp_rat_clear(r);
144
145 free(ptr: r);
146}
147
148mp_result mp_rat_numer(mp_rat r, mp_int z) {
149 return mp_int_copy(a: MP_NUMER_P(Q: r), c: z);
150}
151
152mp_int mp_rat_numer_ref(mp_rat r) { return MP_NUMER_P(Q: r); }
153
154mp_result mp_rat_denom(mp_rat r, mp_int z) {
155 return mp_int_copy(a: MP_DENOM_P(Q: r), c: z);
156}
157
158mp_int mp_rat_denom_ref(mp_rat r) { return MP_DENOM_P(Q: r); }
159
160mp_sign mp_rat_sign(mp_rat r) { return MP_NUMER_SIGN(r); }
161
162mp_result mp_rat_copy(mp_rat a, mp_rat c) {
163 mp_result res;
164
165 if ((res = mp_int_copy(a: MP_NUMER_P(Q: a), c: MP_NUMER_P(Q: c))) != MP_OK) {
166 return res;
167 }
168
169 res = mp_int_copy(a: MP_DENOM_P(Q: a), c: MP_DENOM_P(Q: c));
170 return res;
171}
172
173void mp_rat_zero(mp_rat r) {
174 mp_int_zero(z: MP_NUMER_P(Q: r));
175 mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
176}
177
178mp_result mp_rat_abs(mp_rat a, mp_rat c) {
179 mp_result res;
180
181 if ((res = mp_int_abs(a: MP_NUMER_P(Q: a), c: MP_NUMER_P(Q: c))) != MP_OK) {
182 return res;
183 }
184
185 res = mp_int_abs(a: MP_DENOM_P(Q: a), c: MP_DENOM_P(Q: c));
186 return res;
187}
188
189mp_result mp_rat_neg(mp_rat a, mp_rat c) {
190 mp_result res;
191
192 if ((res = mp_int_neg(a: MP_NUMER_P(Q: a), c: MP_NUMER_P(Q: c))) != MP_OK) {
193 return res;
194 }
195
196 res = mp_int_copy(a: MP_DENOM_P(Q: a), c: MP_DENOM_P(Q: c));
197 return res;
198}
199
200mp_result mp_rat_recip(mp_rat a, mp_rat c) {
201 mp_result res;
202
203 if (mp_rat_compare_zero(r: a) == 0) return MP_UNDEF;
204
205 if ((res = mp_rat_copy(a, c)) != MP_OK) return res;
206
207 mp_int_swap(a: MP_NUMER_P(Q: c), c: MP_DENOM_P(Q: c));
208
209 /* Restore the signs of the swapped elements */
210 {
211 mp_sign tmp = MP_NUMER_SIGN(c);
212
213 MP_NUMER_SIGN(c) = MP_DENOM_SIGN(c);
214 MP_DENOM_SIGN(c) = tmp;
215 }
216
217 return MP_OK;
218}
219
220mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c) {
221 return s_rat_combine(a, b, c, comb_f: mp_int_add);
222}
223
224mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c) {
225 return s_rat_combine(a, b, c, comb_f: mp_int_sub);
226}
227
228mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c) {
229 mp_result res;
230
231 if ((res = mp_int_mul(a: MP_NUMER_P(Q: a), b: MP_NUMER_P(Q: b), c: MP_NUMER_P(Q: c))) != MP_OK)
232 return res;
233
234 if (mp_int_compare_zero(z: MP_NUMER_P(Q: c)) != 0) {
235 res = mp_int_mul(a: MP_DENOM_P(Q: a), b: MP_DENOM_P(Q: b), c: MP_DENOM_P(Q: c));
236 if (res != MP_OK) {
237 return res;
238 }
239 }
240
241 return s_rat_reduce(r: c);
242}
243
244mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c) {
245 mp_result res = MP_OK;
246
247 if (mp_rat_compare_zero(r: b) == 0) return MP_UNDEF;
248
249 if (c == a || c == b) {
250 mpz_t tmp;
251
252 if ((res = mp_int_init(z: &tmp)) != MP_OK) return res;
253 if ((res = mp_int_mul(a: MP_NUMER_P(Q: a), b: MP_DENOM_P(Q: b), c: &tmp)) != MP_OK) {
254 goto CLEANUP;
255 }
256 if ((res = mp_int_mul(a: MP_DENOM_P(Q: a), b: MP_NUMER_P(Q: b), c: MP_DENOM_P(Q: c))) !=
257 MP_OK) {
258 goto CLEANUP;
259 }
260 res = mp_int_copy(a: &tmp, c: MP_NUMER_P(Q: c));
261
262 CLEANUP:
263 mp_int_clear(z: &tmp);
264 } else {
265 if ((res = mp_int_mul(a: MP_NUMER_P(Q: a), b: MP_DENOM_P(Q: b), c: MP_NUMER_P(Q: c))) !=
266 MP_OK) {
267 return res;
268 }
269 if ((res = mp_int_mul(a: MP_DENOM_P(Q: a), b: MP_NUMER_P(Q: b), c: MP_DENOM_P(Q: c))) !=
270 MP_OK) {
271 return res;
272 }
273 }
274
275 if (res != MP_OK) {
276 return res;
277 } else {
278 return s_rat_reduce(r: c);
279 }
280}
281
282mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c) {
283 mpz_t tmp;
284 mp_result res;
285
286 if ((res = mp_int_init_copy(z: &tmp, old: b)) != MP_OK) {
287 return res;
288 }
289 if ((res = mp_int_mul(a: &tmp, b: MP_DENOM_P(Q: a), c: &tmp)) != MP_OK) {
290 goto CLEANUP;
291 }
292 if ((res = mp_rat_copy(a, c)) != MP_OK) {
293 goto CLEANUP;
294 }
295 if ((res = mp_int_add(a: MP_NUMER_P(Q: c), b: &tmp, c: MP_NUMER_P(Q: c))) != MP_OK) {
296 goto CLEANUP;
297 }
298
299 res = s_rat_reduce(r: c);
300
301CLEANUP:
302 mp_int_clear(z: &tmp);
303 return res;
304}
305
306mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c) {
307 mpz_t tmp;
308 mp_result res;
309
310 if ((res = mp_int_init_copy(z: &tmp, old: b)) != MP_OK) {
311 return res;
312 }
313 if ((res = mp_int_mul(a: &tmp, b: MP_DENOM_P(Q: a), c: &tmp)) != MP_OK) {
314 goto CLEANUP;
315 }
316 if ((res = mp_rat_copy(a, c)) != MP_OK) {
317 goto CLEANUP;
318 }
319 if ((res = mp_int_sub(a: MP_NUMER_P(Q: c), b: &tmp, c: MP_NUMER_P(Q: c))) != MP_OK) {
320 goto CLEANUP;
321 }
322
323 res = s_rat_reduce(r: c);
324
325CLEANUP:
326 mp_int_clear(z: &tmp);
327 return res;
328}
329
330mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c) {
331 mp_result res;
332
333 if ((res = mp_rat_copy(a, c)) != MP_OK) {
334 return res;
335 }
336 if ((res = mp_int_mul(a: MP_NUMER_P(Q: c), b, c: MP_NUMER_P(Q: c))) != MP_OK) {
337 return res;
338 }
339
340 return s_rat_reduce(r: c);
341}
342
343mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c) {
344 mp_result res;
345
346 if (mp_int_compare_zero(z: b) == 0) {
347 return MP_UNDEF;
348 }
349 if ((res = mp_rat_copy(a, c)) != MP_OK) {
350 return res;
351 }
352 if ((res = mp_int_mul(a: MP_DENOM_P(Q: c), b, c: MP_DENOM_P(Q: c))) != MP_OK) {
353 return res;
354 }
355
356 return s_rat_reduce(r: c);
357}
358
359mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c) {
360 mp_result res;
361
362 /* Special cases for easy powers. */
363 if (b == 0) {
364 return mp_rat_set_value(r: c, numer: 1, denom: 1);
365 } else if (b == 1) {
366 return mp_rat_copy(a, c);
367 }
368
369 /* Since rationals are always stored in lowest terms, it is not necessary to
370 reduce again when raising to an integer power. */
371 if ((res = mp_int_expt(a: MP_NUMER_P(Q: a), b, c: MP_NUMER_P(Q: c))) != MP_OK) {
372 return res;
373 }
374
375 return mp_int_expt(a: MP_DENOM_P(Q: a), b, c: MP_DENOM_P(Q: c));
376}
377
378int mp_rat_compare(mp_rat a, mp_rat b) {
379 /* Quick check for opposite signs. Works because the sign of the numerator
380 is always definitive. */
381 if (MP_NUMER_SIGN(a) != MP_NUMER_SIGN(b)) {
382 if (MP_NUMER_SIGN(a) == MP_ZPOS) {
383 return 1;
384 } else {
385 return -1;
386 }
387 } else {
388 /* Compare absolute magnitudes; if both are positive, the answer stands,
389 otherwise it needs to be reflected about zero. */
390 int cmp = mp_rat_compare_unsigned(a, b);
391
392 if (MP_NUMER_SIGN(a) == MP_ZPOS) {
393 return cmp;
394 } else {
395 return -cmp;
396 }
397 }
398}
399
400int mp_rat_compare_unsigned(mp_rat a, mp_rat b) {
401 /* If the denominators are equal, we can quickly compare numerators without
402 multiplying. Otherwise, we actually have to do some work. */
403 if (mp_int_compare_unsigned(a: MP_DENOM_P(Q: a), b: MP_DENOM_P(Q: b)) == 0) {
404 return mp_int_compare_unsigned(a: MP_NUMER_P(Q: a), b: MP_NUMER_P(Q: b));
405 }
406
407 else {
408 mpz_t temp[2];
409 mp_result res;
410 int cmp = INT_MAX, last = 0;
411
412 /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
413 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
414 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
415
416 if ((res = mp_int_mul(TEMP(0), b: MP_DENOM_P(Q: b), TEMP(0))) != MP_OK ||
417 (res = mp_int_mul(TEMP(1), b: MP_DENOM_P(Q: a), TEMP(1))) != MP_OK) {
418 goto CLEANUP;
419 }
420
421 cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1));
422
423 CLEANUP:
424 while (--last >= 0) {
425 mp_int_clear(TEMP(last));
426 }
427
428 return cmp;
429 }
430}
431
432int mp_rat_compare_zero(mp_rat r) { return mp_int_compare_zero(z: MP_NUMER_P(Q: r)); }
433
434int mp_rat_compare_value(mp_rat r, mp_small n, mp_small d) {
435 mpq_t tmp;
436 mp_result res;
437 int out = INT_MAX;
438
439 if ((res = mp_rat_init(r: &tmp)) != MP_OK) {
440 return out;
441 }
442 if ((res = mp_rat_set_value(r: &tmp, numer: n, denom: d)) != MP_OK) {
443 goto CLEANUP;
444 }
445
446 out = mp_rat_compare(a: r, b: &tmp);
447
448CLEANUP:
449 mp_rat_clear(r: &tmp);
450 return out;
451}
452
453bool mp_rat_is_integer(mp_rat r) {
454 return (mp_int_compare_value(z: MP_DENOM_P(Q: r), v: 1) == 0);
455}
456
457mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den) {
458 mp_result res;
459
460 if ((res = mp_int_to_int(z: MP_NUMER_P(Q: r), out: num)) != MP_OK) {
461 return res;
462 }
463
464 res = mp_int_to_int(z: MP_DENOM_P(Q: r), out: den);
465 return res;
466}
467
468mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit) {
469 /* Write the numerator. The sign of the rational number is written by the
470 underlying integer implementation. */
471 mp_result res;
472 if ((res = mp_int_to_string(z: MP_NUMER_P(Q: r), radix, str, limit)) != MP_OK) {
473 return res;
474 }
475
476 /* If the value is zero, don't bother writing any denominator */
477 if (mp_int_compare_zero(z: MP_NUMER_P(Q: r)) == 0) {
478 return MP_OK;
479 }
480
481 /* Locate the end of the numerator, and make sure we are not going to exceed
482 the limit by writing a slash. */
483 int len = strlen(s: str);
484 char *start = str + len;
485 limit -= len;
486 if (limit == 0) return MP_TRUNC;
487
488 *start++ = '/';
489 limit -= 1;
490
491 return mp_int_to_string(z: MP_DENOM_P(Q: r), radix, str: start, limit);
492}
493
494mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec,
495 mp_round_mode round, char *str, int limit) {
496 mpz_t temp[3];
497 mp_result res;
498 int last = 0;
499
500 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(r)), last);
501 SETUP(mp_int_init(TEMP(last)), last);
502 SETUP(mp_int_init(TEMP(last)), last);
503
504 /* Get the unsigned integer part by dividing denominator into the absolute
505 value of the numerator. */
506 mp_int_abs(TEMP(0), TEMP(0));
507 if ((res = mp_int_div(TEMP(0), b: MP_DENOM_P(Q: r), TEMP(0), TEMP(1))) != MP_OK) {
508 goto CLEANUP;
509 }
510
511 /* Now: T0 = integer portion, unsigned;
512 T1 = remainder, from which fractional part is computed. */
513
514 /* Count up leading zeroes after the radix point. */
515 int zprec = (int)prec;
516 int lead_0;
517 for (lead_0 = 0; lead_0 < zprec && mp_int_compare(TEMP(1), b: MP_DENOM_P(Q: r)) < 0;
518 ++lead_0) {
519 if ((res = mp_int_mul_value(TEMP(1), value: radix, TEMP(1))) != MP_OK) {
520 goto CLEANUP;
521 }
522 }
523
524 /* Multiply remainder by a power of the radix sufficient to get the right
525 number of significant figures. */
526 if (zprec > lead_0) {
527 if ((res = mp_int_expt_value(a: radix, b: zprec - lead_0, TEMP(2))) != MP_OK) {
528 goto CLEANUP;
529 }
530 if ((res = mp_int_mul(TEMP(1), TEMP(2), TEMP(1))) != MP_OK) {
531 goto CLEANUP;
532 }
533 }
534 if ((res = mp_int_div(TEMP(1), b: MP_DENOM_P(Q: r), TEMP(1), TEMP(2))) != MP_OK) {
535 goto CLEANUP;
536 }
537
538 /* Now: T1 = significant digits of fractional part;
539 T2 = leftovers, to use for rounding.
540
541 At this point, what we do depends on the rounding mode. The default is
542 MP_ROUND_DOWN, for which everything is as it should be already.
543 */
544 switch (round) {
545 int cmp;
546
547 case MP_ROUND_UP:
548 if (mp_int_compare_zero(TEMP(2)) != 0) {
549 if (prec == 0) {
550 res = mp_int_add_value(TEMP(0), value: 1, TEMP(0));
551 } else {
552 res = mp_int_add_value(TEMP(1), value: 1, TEMP(1));
553 }
554 }
555 break;
556
557 case MP_ROUND_HALF_UP:
558 case MP_ROUND_HALF_DOWN:
559 if ((res = mp_int_mul_pow2(TEMP(2), p2: 1, TEMP(2))) != MP_OK) {
560 goto CLEANUP;
561 }
562
563 cmp = mp_int_compare(TEMP(2), b: MP_DENOM_P(Q: r));
564
565 if (round == MP_ROUND_HALF_UP) cmp += 1;
566
567 if (cmp > 0) {
568 if (prec == 0) {
569 res = mp_int_add_value(TEMP(0), value: 1, TEMP(0));
570 } else {
571 res = mp_int_add_value(TEMP(1), value: 1, TEMP(1));
572 }
573 }
574 break;
575
576 case MP_ROUND_DOWN:
577 break; /* No action required */
578
579 default:
580 return MP_BADARG; /* Invalid rounding specifier */
581 }
582 if (res != MP_OK) {
583 goto CLEANUP;
584 }
585
586 /* The sign of the output should be the sign of the numerator, but if all the
587 displayed digits will be zero due to the precision, a negative shouldn't
588 be shown. */
589 char *start = str;
590 int left = limit;
591 if (MP_NUMER_SIGN(r) == MP_NEG && (mp_int_compare_zero(TEMP(0)) != 0 ||
592 mp_int_compare_zero(TEMP(1)) != 0)) {
593 *start++ = '-';
594 left -= 1;
595 }
596
597 if ((res = mp_int_to_string(TEMP(0), radix, str: start, limit: left)) != MP_OK) {
598 goto CLEANUP;
599 }
600
601 int len = strlen(s: start);
602 start += len;
603 left -= len;
604
605 if (prec == 0) goto CLEANUP;
606
607 *start++ = '.';
608 left -= 1;
609
610 if (left < zprec + 1) {
611 res = MP_TRUNC;
612 goto CLEANUP;
613 }
614
615 memset(s: start, c: '0', n: lead_0 - 1);
616 left -= lead_0;
617 start += lead_0 - 1;
618
619 res = mp_int_to_string(TEMP(1), radix, str: start, limit: left);
620
621CLEANUP:
622 while (--last >= 0) mp_int_clear(TEMP(last));
623
624 return res;
625}
626
627mp_result mp_rat_string_len(mp_rat r, mp_size radix) {
628 mp_result d_len = 0;
629 mp_result n_len = mp_int_string_len(z: MP_NUMER_P(Q: r), radix);
630
631 if (mp_int_compare_zero(z: MP_NUMER_P(Q: r)) != 0) {
632 d_len = mp_int_string_len(z: MP_DENOM_P(Q: r), radix);
633 }
634
635 /* Though simplistic, this formula is correct. Space for the sign flag is
636 included in n_len, and the space for the NUL that is counted in n_len
637 counts for the separator here. The space for the NUL counted in d_len
638 counts for the final terminator here. */
639
640 return n_len + d_len;
641}
642
643mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec) {
644 int f_len;
645 int z_len = mp_int_string_len(z: MP_NUMER_P(Q: r), radix);
646
647 if (prec == 0) {
648 f_len = 1; /* terminator only */
649 } else {
650 f_len = 1 + prec + 1; /* decimal point, digits, terminator */
651 }
652
653 return z_len + f_len;
654}
655
656mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str) {
657 return mp_rat_read_cstring(r, radix, str, NULL);
658}
659
660mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str,
661 char **end) {
662 mp_result res;
663 char *endp;
664
665 if ((res = mp_int_read_cstring(z: MP_NUMER_P(Q: r), radix, str, end: &endp)) != MP_OK &&
666 (res != MP_TRUNC)) {
667 return res;
668 }
669
670 /* Skip whitespace between numerator and (possible) separator */
671 while (isspace((unsigned char)*endp)) {
672 ++endp;
673 }
674
675 /* If there is no separator, we will stop reading at this point. */
676 if (*endp != '/') {
677 mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
678 if (end != NULL) *end = endp;
679 return res;
680 }
681
682 ++endp; /* skip separator */
683 if ((res = mp_int_read_cstring(z: MP_DENOM_P(Q: r), radix, str: endp, end)) != MP_OK) {
684 return res;
685 }
686
687 /* Make sure the value is well-defined */
688 if (mp_int_compare_zero(z: MP_DENOM_P(Q: r)) == 0) {
689 return MP_UNDEF;
690 }
691
692 /* Reduce to lowest terms */
693 return s_rat_reduce(r);
694}
695
696/* Read a string and figure out what format it's in. The radix may be supplied
697 as zero to use "default" behaviour.
698
699 This function will accept either a/b notation or decimal notation.
700 */
701mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str,
702 char **end) {
703 char *endp = "";
704 mp_result res;
705
706 if (radix == 0) radix = 10; /* default to decimal input */
707
708 res = mp_rat_read_cstring(r, radix, str, end: &endp);
709 if (res == MP_TRUNC && *endp == '.') {
710 res = mp_rat_read_cdecimal(r, radix, str, end: &endp);
711 }
712 if (end != NULL) *end = endp;
713 return res;
714}
715
716mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str) {
717 return mp_rat_read_cdecimal(r, radix, str, NULL);
718}
719
720mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str,
721 char **end) {
722 mp_result res;
723 mp_sign osign;
724 char *endp;
725
726 while (isspace((unsigned char)*str)) ++str;
727
728 switch (*str) {
729 case '-':
730 osign = MP_NEG;
731 break;
732 default:
733 osign = MP_ZPOS;
734 }
735
736 if ((res = mp_int_read_cstring(z: MP_NUMER_P(Q: r), radix, str, end: &endp)) != MP_OK &&
737 (res != MP_TRUNC)) {
738 return res;
739 }
740
741 /* This needs to be here. */
742 (void)mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
743
744 if (*endp != '.') {
745 if (end != NULL) *end = endp;
746 return res;
747 }
748
749 /* If the character following the decimal point is whitespace or a sign flag,
750 we will consider this a truncated value. This special case is because
751 mp_int_read_string() will consider whitespace or sign flags to be valid
752 starting characters for a value, and we do not want them following the
753 decimal point.
754
755 Once we have done this check, it is safe to read in the value of the
756 fractional piece as a regular old integer.
757 */
758 ++endp;
759 if (*endp == '\0') {
760 if (end != NULL) *end = endp;
761 return MP_OK;
762 } else if (isspace((unsigned char)*endp) || *endp == '-' || *endp == '+') {
763 return MP_TRUNC;
764 } else {
765 mpz_t frac;
766 mp_result save_res;
767 char *save = endp;
768 int num_lz = 0;
769
770 /* Make a temporary to hold the part after the decimal point. */
771 if ((res = mp_int_init(z: &frac)) != MP_OK) {
772 return res;
773 }
774
775 if ((res = mp_int_read_cstring(z: &frac, radix, str: endp, end: &endp)) != MP_OK &&
776 (res != MP_TRUNC)) {
777 goto CLEANUP;
778 }
779
780 /* Save this response for later. */
781 save_res = res;
782
783 if (mp_int_compare_zero(z: &frac) == 0) goto FINISHED;
784
785 /* Discard trailing zeroes (somewhat inefficiently) */
786 while (mp_int_divisible_value(a: &frac, v: radix)) {
787 if ((res = mp_int_div_value(a: &frac, value: radix, q: &frac, NULL)) != MP_OK) {
788 goto CLEANUP;
789 }
790 }
791
792 /* Count leading zeros after the decimal point */
793 while (save[num_lz] == '0') {
794 ++num_lz;
795 }
796
797 /* Find the least power of the radix that is at least as large as the
798 significant value of the fractional part, ignoring leading zeroes. */
799 (void)mp_int_set_value(z: MP_DENOM_P(Q: r), value: radix);
800
801 while (mp_int_compare(a: MP_DENOM_P(Q: r), b: &frac) < 0) {
802 if ((res = mp_int_mul_value(a: MP_DENOM_P(Q: r), value: radix, c: MP_DENOM_P(Q: r))) !=
803 MP_OK) {
804 goto CLEANUP;
805 }
806 }
807
808 /* Also shift by enough to account for leading zeroes */
809 while (num_lz > 0) {
810 if ((res = mp_int_mul_value(a: MP_DENOM_P(Q: r), value: radix, c: MP_DENOM_P(Q: r))) !=
811 MP_OK) {
812 goto CLEANUP;
813 }
814
815 --num_lz;
816 }
817
818 /* Having found this power, shift the numerator leftward that many, digits,
819 and add the nonzero significant digits of the fractional part to get the
820 result. */
821 if ((res = mp_int_mul(a: MP_NUMER_P(Q: r), b: MP_DENOM_P(Q: r), c: MP_NUMER_P(Q: r))) !=
822 MP_OK) {
823 goto CLEANUP;
824 }
825
826 { /* This addition needs to be unsigned. */
827 MP_NUMER_SIGN(r) = MP_ZPOS;
828 if ((res = mp_int_add(a: MP_NUMER_P(Q: r), b: &frac, c: MP_NUMER_P(Q: r))) != MP_OK) {
829 goto CLEANUP;
830 }
831
832 MP_NUMER_SIGN(r) = osign;
833 }
834 if ((res = s_rat_reduce(r)) != MP_OK) goto CLEANUP;
835
836 /* At this point, what we return depends on whether reading the fractional
837 part was truncated or not. That information is saved from when we
838 called mp_int_read_string() above. */
839 FINISHED:
840 res = save_res;
841 if (end != NULL) *end = endp;
842
843 CLEANUP:
844 mp_int_clear(z: &frac);
845
846 return res;
847 }
848}
849
850/* Private functions for internal use. Make unchecked assumptions about format
851 and validity of inputs. */
852
853static mp_result s_rat_reduce(mp_rat r) {
854 mpz_t gcd;
855 mp_result res = MP_OK;
856
857 if (mp_int_compare_zero(z: MP_NUMER_P(Q: r)) == 0) {
858 mp_int_set_value(z: MP_DENOM_P(Q: r), value: 1);
859 return MP_OK;
860 }
861
862 /* If the greatest common divisor of the numerator and denominator is greater
863 than 1, divide it out. */
864 if ((res = mp_int_init(z: &gcd)) != MP_OK) return res;
865
866 if ((res = mp_int_gcd(a: MP_NUMER_P(Q: r), b: MP_DENOM_P(Q: r), c: &gcd)) != MP_OK) {
867 goto CLEANUP;
868 }
869
870 if (mp_int_compare_value(z: &gcd, v: 1) != 0) {
871 if ((res = mp_int_div(a: MP_NUMER_P(Q: r), b: &gcd, q: MP_NUMER_P(Q: r), NULL)) != MP_OK) {
872 goto CLEANUP;
873 }
874 if ((res = mp_int_div(a: MP_DENOM_P(Q: r), b: &gcd, q: MP_DENOM_P(Q: r), NULL)) != MP_OK) {
875 goto CLEANUP;
876 }
877 }
878
879 /* Fix up the signs of numerator and denominator */
880 if (MP_NUMER_SIGN(r) == MP_DENOM_SIGN(r)) {
881 MP_NUMER_SIGN(r) = MP_DENOM_SIGN(r) = MP_ZPOS;
882 } else {
883 MP_NUMER_SIGN(r) = MP_NEG;
884 MP_DENOM_SIGN(r) = MP_ZPOS;
885 }
886
887CLEANUP:
888 mp_int_clear(z: &gcd);
889
890 return res;
891}
892
893static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,
894 mp_result (*comb_f)(mp_int, mp_int, mp_int)) {
895 mp_result res;
896
897 /* Shortcut when denominators are already common */
898 if (mp_int_compare(a: MP_DENOM_P(Q: a), b: MP_DENOM_P(Q: b)) == 0) {
899 if ((res = (comb_f)(MP_NUMER_P(Q: a), MP_NUMER_P(Q: b), MP_NUMER_P(Q: c))) !=
900 MP_OK) {
901 return res;
902 }
903 if ((res = mp_int_copy(a: MP_DENOM_P(Q: a), c: MP_DENOM_P(Q: c))) != MP_OK) {
904 return res;
905 }
906
907 return s_rat_reduce(r: c);
908 } else {
909 mpz_t temp[2];
910 int last = 0;
911
912 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
913 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
914
915 if ((res = mp_int_mul(TEMP(0), b: MP_DENOM_P(Q: b), TEMP(0))) != MP_OK) {
916 goto CLEANUP;
917 }
918 if ((res = mp_int_mul(TEMP(1), b: MP_DENOM_P(Q: a), TEMP(1))) != MP_OK) {
919 goto CLEANUP;
920 }
921 if ((res = (comb_f)(TEMP(0), TEMP(1), MP_NUMER_P(Q: c))) != MP_OK) {
922 goto CLEANUP;
923 }
924
925 res = mp_int_mul(a: MP_DENOM_P(Q: a), b: MP_DENOM_P(Q: b), c: MP_DENOM_P(Q: c));
926
927 CLEANUP:
928 while (--last >= 0) {
929 mp_int_clear(TEMP(last));
930 }
931
932 if (res == MP_OK) {
933 return s_rat_reduce(r: c);
934 } else {
935 return res;
936 }
937 }
938}
939
940/* Here there be dragons */
941

source code of polly/lib/External/isl/imath/imrat.c