# klisp

an open source interpreter for the Kernel Programming Language.
git clone http://git.hanabi.in/repos/klisp.git
Log | Files | Refs | README

imrat.c (27782B)

```      1 /*
2 ** imrat.c
3 ** Arbitrary precision rational arithmetic routines.
4 ** See Copyright Notice in klisp.h
5 */
6
7 /*
8 ** SOURCE NOTE: This is mostly from the IMath library, written by
9 ** M.J. Fromberger. It is adapted to klisp, mainly in the use of the
10 ** klisp allocator and fixing of digit size to 32 bits.
11 ** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM.
12 */
13
14 #include "imrat.h"
15 #include <stdlib.h>
16 #include <string.h>
17 #include <ctype.h>
18 #include <assert.h>
19
20 /* Andres Navarro: klisp includes */
21 #include "kobject.h"
22 #include "kstate.h"
23 #include "kmem.h"
24 #include "kerror.h"
25
26 /* {{{ Useful macros */
27
28 #define TEMP(K) (temp + (K))
29 #define SETUP(E, C)                                             \
30     do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
31
32 /* Argument checking:
33    Use CHECK() where a return value is required; NRCHECK() elsewhere */
34 #define CHECK(TEST)   assert(TEST)
35 #define NRCHECK(TEST) assert(TEST)
36
37 /* }}} */
38
39 /* Reduce the given rational, in place, to lowest terms and canonical
40    form.  Zero is represented as 0/1, one as 1/1.  Signs are adjusted
41    so that the sign of the numerator is definitive. */
42 static mp_result s_rat_reduce(klisp_State *K, mp_rat r);
43
44 /* Common code for addition and subtraction operations on rationals. */
45 static mp_result s_rat_combine(klisp_State *K, mp_rat a, mp_rat b, mp_rat c,
46                                mp_result (*comb_f)
47                                (klisp_State *,mp_int, mp_int, mp_int));
48
49 /* {{{ mp_rat_init(r) */
50
51 mp_result mp_rat_init(klisp_State *K, mp_rat r)
52 {
53     return mp_rat_init_size(K, r, 0, 0);
54 }
55
56 /* }}} */
57
58 /* {{{ mp_rat_alloc() */
59
60 mp_rat mp_rat_alloc(klisp_State *K)
61 {
62     mp_rat out = klispM_new(K, mpq_t);
63
64     if(out != NULL) {
65         if(mp_rat_init(K, out) != MP_OK) {
66             klispM_free(K, out);
67             return NULL;
68         }
69     }
70
71     return out;
72 }
73
74 /* }}} */
75
76 /* {{{ mp_rat_init_size(r, n_prec, d_prec) */
77
78 mp_result mp_rat_init_size(klisp_State *K, mp_rat r, mp_size n_prec,
79                            mp_size d_prec)
80 {
81     mp_result res;
82
83     if((res = mp_int_init_size(K, MP_NUMER_P(r), n_prec)) != MP_OK)
84         return res;
85     if((res = mp_int_init_size(K, MP_DENOM_P(r), d_prec)) != MP_OK) {
86         mp_int_clear(K, MP_NUMER_P(r));
87         return res;
88     }
89
90     return mp_int_set_value(K, MP_DENOM_P(r), 1);
91 }
92
93 /* }}} */
94
95 /* {{{ mp_rat_init_copy(r, old) */
96
97 mp_result mp_rat_init_copy(klisp_State *K, mp_rat r, mp_rat old)
98 {
99     mp_result res;
100
101     if((res = mp_int_init_copy(K, MP_NUMER_P(r), MP_NUMER_P(old))) != MP_OK)
102         return res;
103     if((res = mp_int_init_copy(K, MP_DENOM_P(r), MP_DENOM_P(old))) != MP_OK)
104         mp_int_clear(K, MP_NUMER_P(r));
105
106     return res;
107 }
108
109 /* }}} */
110
111 /* {{{ mp_rat_set_value(r, numer, denom) */
112
113 mp_result mp_rat_set_value(klisp_State *K, mp_rat r, int numer, int denom)
114 {
115     mp_result res;
116
117     if(denom == 0)
118         return MP_UNDEF;
119
120     if((res = mp_int_set_value(K, MP_NUMER_P(r), numer)) != MP_OK)
121         return res;
122     if((res = mp_int_set_value(K, MP_DENOM_P(r), denom)) != MP_OK)
123         return res;
124
125     return s_rat_reduce(K, r);
126 }
127
128 /* }}} */
129
130 /* {{{ mp_rat_clear(r) */
131
132 void         mp_rat_clear(klisp_State *K, mp_rat r)
133 {
134     mp_int_clear(K, MP_NUMER_P(r));
135     mp_int_clear(K, MP_DENOM_P(r));
136 }
137
138 /* }}} */
139
140 /* {{{ mp_rat_free(r) */
141
142 void         mp_rat_free(klisp_State *K, mp_rat r)
143 {
144     NRCHECK(r != NULL);
145
146     if(r->num.digits != NULL)
147         mp_rat_clear(K, r);
148
149     klispM_free(K, r);
150 }
151
152 /* }}} */
153
154 /* {{{ mp_rat_numer(r, z) */
155
156 mp_result mp_rat_numer(klisp_State *K, mp_rat r, mp_int z)
157 {
158     return mp_int_copy(K, MP_NUMER_P(r), z);
159 }
160
161 /* }}} */
162
163 /* {{{ mp_rat_denom(r, z) */
164
165 mp_result mp_rat_denom(klisp_State *K, mp_rat r, mp_int z)
166 {
167     return mp_int_copy(K, MP_DENOM_P(r), z);
168 }
169
170 /* }}} */
171
172 /* {{{ mp_rat_sign(r) */
173
174 mp_sign   mp_rat_sign(mp_rat r)
175 {
176     return MP_SIGN(MP_NUMER_P(r));
177 }
178
179 /* }}} */
180
181 /* {{{ mp_rat_copy(a, c) */
182
183 mp_result mp_rat_copy(klisp_State *K, mp_rat a, mp_rat c)
184 {
185     mp_result res;
186
187     if((res = mp_int_copy(K, MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
188         return res;
189
190     res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c));
191     return res;
192 }
193
194 /* }}} */
195
196 /* {{{ mp_rat_zero(r) */
197
198 void         mp_rat_zero(klisp_State *K, mp_rat r)
199 {
200     mp_int_zero(MP_NUMER_P(r));
201     mp_int_set_value(K, MP_DENOM_P(r), 1);
202 }
203
204 /* }}} */
205
206 /* {{{ mp_rat_abs(a, c) */
207
208 mp_result mp_rat_abs(klisp_State *K, mp_rat a, mp_rat c)
209 {
210     mp_result res;
211
212     if((res = mp_int_abs(K, MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
213         return res;
214
215     res = mp_int_abs(K, MP_DENOM_P(a), MP_DENOM_P(c));
216     return res;
217 }
218
219 /* }}} */
220
221 /* {{{ mp_rat_neg(a, c) */
222
223 mp_result mp_rat_neg(klisp_State *K, mp_rat a, mp_rat c)
224 {
225     mp_result res;
226
227     if((res = mp_int_neg(K, MP_NUMER_P(a),
228                          MP_NUMER_P(c))) != MP_OK)
229         return res;
230
231     res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c));
232     return res;
233 }
234
235 /* }}} */
236
237 /* {{{ mp_rat_recip(a, c) */
238
239 mp_result mp_rat_recip(klisp_State *K, mp_rat a, mp_rat c)
240 {
241     mp_result res;
242
243     if(mp_rat_compare_zero(a) == 0)
244         return MP_UNDEF;
245
246     if((res = mp_rat_copy(K, a, c)) != MP_OK)
247         return res;
248
249     mp_int_swap(MP_NUMER_P(c), MP_DENOM_P(c));
250
251     /* Restore the signs of the swapped elements */
252     {
253         mp_sign tmp = MP_SIGN(MP_NUMER_P(c));
254
255         MP_SIGN(MP_NUMER_P(c)) = MP_SIGN(MP_DENOM_P(c));
256         MP_SIGN(MP_DENOM_P(c)) = tmp;
257     }
258
259     return MP_OK;
260 }
261
262 /* }}} */
263
264 /* {{{ mp_rat_add(a, b, c) */
265
266 mp_result mp_rat_add(klisp_State *K, mp_rat a, mp_rat b, mp_rat c)
267 {
268     return s_rat_combine(K, a, b, c, mp_int_add);
269
270 }
271
272 /* }}} */
273
274 /* {{{ mp_rat_sub(a, b, c) */
275
276 mp_result mp_rat_sub(klisp_State *K, mp_rat a, mp_rat b, mp_rat c)
277 {
278     return s_rat_combine(K, a, b, c, mp_int_sub);
279
280 }
281
282 /* }}} */
283
284 /* {{{ mp_rat_mul(a, b, c) */
285
286 mp_result mp_rat_mul(klisp_State *K, mp_rat a, mp_rat b, mp_rat c)
287 {
288     mp_result res;
289
290     if((res = mp_int_mul(K, MP_NUMER_P(a), MP_NUMER_P(b),
291                          MP_NUMER_P(c))) != MP_OK)
292         return res;
293
294     if(mp_int_compare_zero(MP_NUMER_P(c)) != 0) {
295         if((res = mp_int_mul(K, MP_DENOM_P(a), MP_DENOM_P(b),
296                              MP_DENOM_P(c))) != MP_OK)
297             return res;
298     }
299
300     return s_rat_reduce(K, c);
301 }
302
303 /* }}} */
304
305 /* {{{ mp_int_div(a, b, c) */
306
307 mp_result mp_rat_div(klisp_State *K, mp_rat a, mp_rat b, mp_rat c)
308 {
309     mp_result res = MP_OK;
310
311     if(mp_rat_compare_zero(b) == 0)
312         return MP_UNDEF;
313
314     if(c == a || c == b) {
315         mpz_t tmp;
316
317         if((res = mp_int_init(&tmp)) != MP_OK) return res;
318         if((res = mp_int_mul(K, MP_NUMER_P(a), MP_DENOM_P(b), &tmp)) != MP_OK)
319             goto CLEANUP;
320         if((res = mp_int_mul(K, MP_DENOM_P(a), MP_NUMER_P(b),
321                              MP_DENOM_P(c))) != MP_OK)
322             goto CLEANUP;
323         res = mp_int_copy(K, &tmp, MP_NUMER_P(c));
324
325     CLEANUP:
326         mp_int_clear(K, &tmp);
327     }
328     else {
329         if((res = mp_int_mul(K, MP_NUMER_P(a), MP_DENOM_P(b),
330                              MP_NUMER_P(c))) != MP_OK)
331             return res;
332         if((res = mp_int_mul(K, MP_DENOM_P(a), MP_NUMER_P(b),
333                              MP_DENOM_P(c))) != MP_OK)
334             return res;
335     }
336
337     if(res != MP_OK)
338         return res;
339     else
340         return s_rat_reduce(K, c);
341 }
342
343 /* }}} */
344
345 /* {{{ mp_rat_add_int(a, b, c) */
346
347 mp_result mp_rat_add_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c)
348 {
349     mpz_t tmp;
350     mp_result res;
351
352     if((res = mp_int_init_copy(K, &tmp, b)) != MP_OK)
353         return res;
354
355     if((res = mp_int_mul(K, &tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
356         goto CLEANUP;
357
358     if((res = mp_rat_copy(K, a, c)) != MP_OK)
359         goto CLEANUP;
360
361     if((res = mp_int_add(K, MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
362         goto CLEANUP;
363
364     res = s_rat_reduce(K, c);
365
366 CLEANUP:
367     mp_int_clear(K, &tmp);
368     return res;
369 }
370
371 /* }}} */
372
373 /* {{{ mp_rat_sub_int(a, b, c) */
374
375 mp_result mp_rat_sub_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c)
376 {
377     mpz_t tmp;
378     mp_result res;
379
380     if((res = mp_int_init_copy(K, &tmp, b)) != MP_OK)
381         return res;
382
383     if((res = mp_int_mul(K, &tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
384         goto CLEANUP;
385
386     if((res = mp_rat_copy(K, a, c)) != MP_OK)
387         goto CLEANUP;
388
389     if((res = mp_int_sub(K, MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
390         goto CLEANUP;
391
392     res = s_rat_reduce(K, c);
393
394 CLEANUP:
395     mp_int_clear(K, &tmp);
396     return res;
397 }
398
399 /* }}} */
400
401 /* {{{ mp_rat_mul_int(a, b, c) */
402
403 mp_result mp_rat_mul_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c)
404 {
405     mp_result res;
406
407     if((res = mp_rat_copy(K, a, c)) != MP_OK)
408         return res;
409
410     if((res = mp_int_mul(K, MP_NUMER_P(c), b, MP_NUMER_P(c))) != MP_OK)
411         return res;
412
413     return s_rat_reduce(K, c);
414 }
415
416 /* }}} */
417
418 /* {{{ mp_rat_div_int(a, b, c) */
419
420 mp_result mp_rat_div_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c)
421 {
422     mp_result res;
423
424     if(mp_int_compare_zero(b) == 0)
425         return MP_UNDEF;
426
427     if((res = mp_rat_copy(K, a, c)) != MP_OK)
428         return res;
429
430     if((res = mp_int_mul(K, MP_DENOM_P(c), b, MP_DENOM_P(c))) != MP_OK)
431         return res;
432
433     return s_rat_reduce(K, c);
434 }
435
436 /* }}} */
437
438 /* {{{ mp_rat_expt(a, b, c) */
439
440 mp_result mp_rat_expt(klisp_State *K, mp_rat a, mp_small b, mp_rat c)
441 {
442     mp_result  res;
443
444     /* Special cases for easy powers. */
445     if(b == 0)
446         return mp_rat_set_value(K, c, 1, 1);
447     else if(b == 1)
448         return mp_rat_copy(K, a, c);
449
450     /* Since rationals are always stored in lowest terms, it is not
451        necessary to reduce again when raising to an integer power. */
452     if((res = mp_int_expt(K, MP_NUMER_P(a), b, MP_NUMER_P(c))) != MP_OK)
453         return res;
454
455     return mp_int_expt(K, MP_DENOM_P(a), b, MP_DENOM_P(c));
456 }
457
458 /* }}} */
459
460 /* {{{ mp_rat_compare(a, b) */
461
462 int          mp_rat_compare(klisp_State *K, mp_rat a, mp_rat b)
463 {
464     /* Quick check for opposite signs.  Works because the sign of the
465        numerator is always definitive. */
466     if(MP_SIGN(MP_NUMER_P(a)) != MP_SIGN(MP_NUMER_P(b))) {
467         if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
468             return 1;
469         else
470             return -1;
471     }
472     else {
473         /* Compare absolute magnitudes; if both are positive, the answer
474            stands, otherwise it needs to be reflected about zero. */
475         int cmp = mp_rat_compare_unsigned(K, a, b);
476
477         if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
478             return cmp;
479         else
480             return -cmp;
481     }
482 }
483
484 /* }}} */
485
486 /* {{{ mp_rat_compare_unsigned(a, b) */
487
488 int          mp_rat_compare_unsigned(klisp_State *K, mp_rat a, mp_rat b)
489 {
490     /* If the denominators are equal, we can quickly compare numerators
491        without multiplying.  Otherwise, we actually have to do some work. */
492     if(mp_int_compare_unsigned(MP_DENOM_P(a), MP_DENOM_P(b)) == 0)
493         return mp_int_compare_unsigned(MP_NUMER_P(a), MP_NUMER_P(b));
494
495     else {
496         mpz_t  temp[2];
497         mp_result res;
498         int  cmp = INT_MAX, last = 0;
499
500         /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
501         SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(a)), last);
502         SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(b)), last);
503
504         if((res = mp_int_mul(K, TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK ||
505            (res = mp_int_mul(K, TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
506             goto CLEANUP;
507
508         cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1));
509
510     CLEANUP:
511         while(--last >= 0)
512             mp_int_clear(K, TEMP(last));
513
514         return cmp;
515     }
516 }
517
518 /* }}} */
519
520 /* {{{ mp_rat_compare_zero(r) */
521
522 int          mp_rat_compare_zero(mp_rat r)
523 {
524     return mp_int_compare_zero(MP_NUMER_P(r));
525
526 }
527
528 /* }}} */
529
530 /* {{{ mp_rat_compare_value(r, n, d) */
531
532 int mp_rat_compare_value(klisp_State *K, mp_rat r, mp_small n, mp_small d)
533 {
534     mpq_t tmp;
535     mp_result res;
536     int  out = INT_MAX;
537
538     if((res = mp_rat_init(K, &tmp)) != MP_OK)
539         return out;
540     if((res = mp_rat_set_value(K, &tmp, n, d)) != MP_OK)
541         goto CLEANUP;
542
543     out = mp_rat_compare(K, r, &tmp);
544
545 CLEANUP:
546     mp_rat_clear(K, &tmp);
547     return out;
548 }
549
550 /* }}} */
551
552 /* {{{ mp_rat_is_integer(r) */
553
554 int          mp_rat_is_integer(mp_rat r)
555 {
556     return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0);
557 }
558
559 /* }}} */
560
561 /* {{{ mp_rat_to_ints(r, *num, *den) */
562
563 mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den)
564 {
565     mp_result res;
566
567     if((res = mp_int_to_int(MP_NUMER_P(r), num)) != MP_OK)
568         return res;
569
570     res = mp_int_to_int(MP_DENOM_P(r), den);
571     return res;
572 }
573
574 /* }}} */
575
576 /* {{{ mp_rat_to_string(r, radix, *str, limit) */
577
578 mp_result mp_rat_to_string(klisp_State *K, mp_rat r, mp_size radix, char *str,
579                            int limit)
580 {
581     char *start;
582     int   len;
583     mp_result res;
584
585     /* Write the numerator.  The sign of the rational number is written
586        by the underlying integer implementation. */
587     if((res = mp_int_to_string(K, MP_NUMER_P(r), radix, str, limit)) != MP_OK)
588         return res;
589
590     /* If the value is zero, don't bother writing any denominator */
591     if(mp_int_compare_zero(MP_NUMER_P(r)) == 0)
592         return MP_OK;
593
594     /* Locate the end of the numerator, and make sure we are not going to
595        exceed the limit by writing a slash. */
596     len = strlen(str);
597     start = str + len;
598     limit -= len;
599     if(limit == 0)
600         return MP_TRUNC;
601
602     *start++ = '/';
603     limit -= 1;
604
605     res = mp_int_to_string(K, MP_DENOM_P(r), radix, start, limit);
606     return res;
607 }
608
609 /* }}} */
610
611 /* {{{ mp_rat_to_decimal(r, radix, prec, *str, limit) */
612 mp_result mp_rat_to_decimal(klisp_State *K, mp_rat r, mp_size radix,
613                             mp_size prec, mp_round_mode round, char *str,
614                             int limit)
615 {
616     mpz_t temp[3];
617     mp_result res;
618     char *start = str;
619     int len, lead_0, left = limit, last = 0;
620
621     SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(r)), last);
622     SETUP(mp_int_init(TEMP(last)), last);
623     SETUP(mp_int_init(TEMP(last)), last);
624
625     /* Get the unsigned integer part by dividing denominator into the
626        absolute value of the numerator. */
627     mp_int_abs(K, TEMP(0), TEMP(0));
628     if((res = mp_int_div(K, TEMP(0), MP_DENOM_P(r), TEMP(0),
629                          TEMP(1))) != MP_OK)
630         goto CLEANUP;
631
632     /* Now:  T0 = integer portion, unsigned;
633        T1 = remainder, from which fractional part is computed. */
634
636     for(lead_0 = 0; lead_0 < prec && mp_int_compare(TEMP(1), MP_DENOM_P(r)) < 0;
638         if((res = mp_int_mul_value(K, TEMP(1), radix, TEMP(1))) != MP_OK)
639             goto CLEANUP;
640     }
641
642     /* Multiply remainder by a power of the radix sufficient to get the
643        right number of significant figures. */
646                                     TEMP(2))) != MP_OK)
647             goto CLEANUP;
648         if((res = mp_int_mul(K, TEMP(1), TEMP(2), TEMP(1))) != MP_OK)
649             goto CLEANUP;
650     }
651     if((res = mp_int_div(K, TEMP(1), MP_DENOM_P(r), TEMP(1),
652                          TEMP(2))) != MP_OK)
653         goto CLEANUP;
654
655     /* Now:  T1 = significant digits of fractional part;
656        T2 = leftovers, to use for rounding.
657
658        At this point, what we do depends on the rounding mode.  The
659        default is MP_ROUND_DOWN, for which everything is as it should be
661     */
662     switch(round) {
663         int cmp;
664
665     case MP_ROUND_UP:
666         if(mp_int_compare_zero(TEMP(2)) != 0) {
667             if(prec == 0)
668                 res = mp_int_add_value(K, TEMP(0), 1, TEMP(0));
669             else
670                 res = mp_int_add_value(K, TEMP(1), 1, TEMP(1));
671         }
672         break;
673
674     case MP_ROUND_HALF_UP:
675     case MP_ROUND_HALF_DOWN:
676         if((res = mp_int_mul_pow2(K, TEMP(2), 1, TEMP(2))) != MP_OK)
677             goto CLEANUP;
678
679         cmp = mp_int_compare(TEMP(2), MP_DENOM_P(r));
680
681         if(round == MP_ROUND_HALF_UP)
682             cmp += 1;
683
684         if(cmp > 0) {
685             if(prec == 0)
686                 res = mp_int_add_value(K, TEMP(0), 1, TEMP(0));
687             else
688                 res = mp_int_add_value(K, TEMP(1), 1, TEMP(1));
689         }
690         break;
691
692     case MP_ROUND_DOWN:
693         break;  /* No action required */
694
695     default:
696         return MP_BADARG; /* Invalid rounding specifier */
697     }
698
699     /* The sign of the output should be the sign of the numerator, but
700        if all the displayed digits will be zero due to the precision, a
701        negative shouldn't be shown. */
702     if(MP_SIGN(MP_NUMER_P(r)) == MP_NEG &&
703        (mp_int_compare_zero(TEMP(0)) != 0 ||
704         mp_int_compare_zero(TEMP(1)) != 0)) {
705         *start++ = '-';
706         left -= 1;
707     }
708
709     if((res = mp_int_to_string(K, TEMP(0), radix, start, left)) != MP_OK)
710         goto CLEANUP;
711
712     len = strlen(start);
713     start += len;
714     left -= len;
715
716     if(prec == 0)
717         goto CLEANUP;
718
719     *start++ = '.';
720     left -= 1;
721
722     if(left < prec + 1) {
723         res = MP_TRUNC;
724         goto CLEANUP;
725     }
726
727     memset(start, '0', lead_0 - 1);
729     start += lead_0 - 1;
730
731     res = mp_int_to_string(K, TEMP(1), radix, start, left);
732
733 CLEANUP:
734     while(--last >= 0)
735         mp_int_clear(K, TEMP(last));
736
737     return res;
738 }
739
740 /* }}} */
741
742 /* {{{ mp_rat_string_len(r, radix) */
743
744 mp_result mp_rat_string_len(mp_rat r, mp_size radix)
745 {
746     mp_result n_len, d_len = 0;
747
749
750     if(mp_int_compare_zero(MP_NUMER_P(r)) != 0)
752
753     /* Though simplistic, this formula is correct.  Space for the sign
754        flag is included in n_len, and the space for the NUL that is
755        counted in n_len counts for the separator here.  The space for
756        the NUL counted in d_len counts for the final terminator here. */
757
758     return n_len + d_len;
759
760 }
761
762 /* }}} */
763
764 /* {{{ mp_rat_decimal_len(r, radix, prec) */
765
766 mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec)
767 {
768     int  z_len, f_len;
769
771
772     if(prec == 0)
773         f_len = 1; /* terminator only */
774     else
775         f_len = 1 + prec + 1; /* decimal point, digits, terminator */
776
777     return z_len + f_len;
778 }
779
780 /* }}} */
781
783
785                              const char *str)
786 {
788 }
789
790 /* }}} */
791
793
795                               const char *str, char **end)
796 {
797     mp_result res;
798     char *endp;
799
801                                   &endp)) != MP_OK && (res != MP_TRUNC))
802         return res;
803
804     /* Skip whitespace between numerator and (possible) separator */
805     while(isspace((unsigned char) *endp))
806         ++endp;
807
808     /* If there is no separator, we will stop reading at this point. */
809     if(*endp != '/') {
810         mp_int_set_value(K, MP_DENOM_P(r), 1);
811         if(end != NULL)
812             *end = endp;
813         return res;
814     }
815
816     ++endp; /* skip separator */
818         return res;
819
820     /* Make sure the value is well-defined */
821     if(mp_int_compare_zero(MP_DENOM_P(r)) == 0)
822         return MP_UNDEF;
823
824     /* Reduce to lowest terms */
825     return s_rat_reduce(K, r);
826 }
827
828 /* }}} */
829
831
832 /* Read a string and figure out what format it's in.  The radix may be
833    supplied as zero to use "default" behaviour.
834
835    This function will accept either a/b notation or decimal notation.
836 */
838                               const char *str, char **end)
839 {
840     char         *endp;
841     mp_result  res;
842
844         radix = 10;  /* default to decimal input */
845
847         if(res == MP_TRUNC) {
848             if(*endp == '.')
850         }
851         else
852             return res;
853     }
854
855     if(end != NULL)
856         *end = endp;
857
858     return res;
859 }
860
861 /* }}} */
862
864
866                               const char *str)
867 {
869 }
870
871 /* }}} */
872
874
876                                const char *str, char **end)
877 {
878     mp_result res;
879     mp_sign   osign;
880     char *endp;
881
882     while(isspace((unsigned char) *str))
883         ++str;
884
885     switch(*str) {
886     case '-':
887         osign = MP_NEG;
888         break;
889     default:
890         osign = MP_ZPOS;
891     }
892
894                                   &endp)) != MP_OK && (res != MP_TRUNC))
895         return res;
896
897     /* This needs to be here. */
898     (void) mp_int_set_value(K, MP_DENOM_P(r), 1);
899
900     if(*endp != '.') {
901         if(end != NULL)
902             *end = endp;
903         return res;
904     }
905
906     /* If the character following the decimal point is whitespace or a
907        sign flag, we will consider this a truncated value.  This special
908        case is because mp_int_read_string() will consider whitespace or
909        sign flags to be valid starting characters for a value, and we do
910        not want them following the decimal point.
911
912        Once we have done this check, it is safe to read in the value of
913        the fractional piece as a regular old integer.
914     */
915     ++endp;
916     if(*endp == '\0') {
917         if(end != NULL)
918             *end = endp;
919         return MP_OK;
920     }
921     else if(isspace((unsigned char) *endp) || *endp == '-' || *endp == '+') {
922         return MP_TRUNC;
923     }
924     else {
925         mpz_t  frac;
926         mp_result save_res;
927         char  *save = endp;
928         int    num_lz = 0;
929
930         /* Make a temporary to hold the part after the decimal point. */
931         if((res = mp_int_init(&frac)) != MP_OK)
932             return res;
933
935            (res != MP_TRUNC))
936             goto CLEANUP;
937
938         /* Save this response for later. */
939         save_res = res;
940
941         if(mp_int_compare_zero(&frac) == 0)
942             goto FINISHED;
943
944         /* Discard trailing zeroes (somewhat inefficiently) */
946             if((res = mp_int_div_value(K, &frac, radix, &frac, NULL)) != MP_OK)
947                 goto CLEANUP;
948
949         /* Count leading zeros after the decimal point */
950         while(save[num_lz] == '0')
951             ++num_lz;
952
953         /* Find the least power of the radix that is at least as large as
954            the significant value of the fractional part, ignoring leading
955            zeroes.  */
957
958         while(mp_int_compare(MP_DENOM_P(r), &frac) < 0) {
959             if((res = mp_int_mul_value(K, MP_DENOM_P(r), radix,
960                                        MP_DENOM_P(r))) != MP_OK)
961                 goto CLEANUP;
962         }
963
964         /* Also shift by enough to account for leading zeroes */
965         while(num_lz > 0) {
966             if((res = mp_int_mul_value(K, MP_DENOM_P(r), radix,
967                                        MP_DENOM_P(r))) != MP_OK)
968                 goto CLEANUP;
969
970             --num_lz;
971         }
972
973         /* Having found this power, shift the numerator leftward that
974            many, digits, and add the nonzero significant digits of the
975            fractional part to get the result. */
976         if((res = mp_int_mul(K, MP_NUMER_P(r), MP_DENOM_P(r),
977                              MP_NUMER_P(r))) != MP_OK)
978             goto CLEANUP;
979
980         { /* This addition needs to be unsigned. */
981             MP_SIGN(MP_NUMER_P(r)) = MP_ZPOS;
982             if((res = mp_int_add(K, MP_NUMER_P(r), &frac,
983                                  MP_NUMER_P(r))) != MP_OK)
984                 goto CLEANUP;
985
986             MP_SIGN(MP_NUMER_P(r)) = osign;
987         }
988         if((res = s_rat_reduce(K, r)) != MP_OK)
989             goto CLEANUP;
990
991         /* At this point, what we return depends on whether reading the
992            fractional part was truncated or not.  That information is
993            saved from when we called mp_int_read_string() above. */
994     FINISHED:
995         res = save_res;
996         if(end != NULL)
997             *end = endp;
998
999     CLEANUP:
1000         mp_int_clear(K, &frac);
1001
1002         return res;
1003     }
1004 }
1005
1006 /* }}} */
1007
1008 /* Private functions for internal use.  Make unchecked assumptions
1009    about format and validity of inputs. */
1010
1011 /* {{{ s_rat_reduce(r) */
1012
1013 static mp_result s_rat_reduce(klisp_State *K, mp_rat r)
1014 {
1015     mpz_t gcd;
1016     mp_result res = MP_OK;
1017
1018     if(mp_int_compare_zero(MP_NUMER_P(r)) == 0) {
1019         mp_int_set_value(K, MP_DENOM_P(r), 1);
1020         return MP_OK;
1021     }
1022
1023     /* If the greatest common divisor of the numerator and denominator
1024        is greater than 1, divide it out. */
1025     if((res = mp_int_init(&gcd)) != MP_OK)
1026         return res;
1027
1028     if((res = mp_int_gcd(K, MP_NUMER_P(r), MP_DENOM_P(r), &gcd)) != MP_OK)
1029         goto CLEANUP;
1030
1031     if(mp_int_compare_value(&gcd, 1) != 0) {
1032         if((res = mp_int_div(K, MP_NUMER_P(r), &gcd, MP_NUMER_P(r),
1033                              NULL)) != MP_OK)
1034             goto CLEANUP;
1035         if((res = mp_int_div(K, MP_DENOM_P(r), &gcd, MP_DENOM_P(r),
1036                              NULL)) != MP_OK)
1037             goto CLEANUP;
1038     }
1039
1040     /* Fix up the signs of numerator and denominator */
1041     if(MP_SIGN(MP_NUMER_P(r)) == MP_SIGN(MP_DENOM_P(r)))
1042         MP_SIGN(MP_NUMER_P(r)) = MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
1043     else {
1044         MP_SIGN(MP_NUMER_P(r)) = MP_NEG;
1045         MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
1046     }
1047
1048 CLEANUP:
1049     mp_int_clear(K, &gcd);
1050
1051     return res;
1052 }
1053
1054 /* }}} */
1055
1056 /* {{{ s_rat_combine(a, b, c, comb_f) */
1057
1058 static mp_result s_rat_combine(klisp_State *K, mp_rat a, mp_rat b, mp_rat c,
1059                                mp_result (*comb_f)(klisp_State *K, mp_int,
1060                                                    mp_int, mp_int))
1061 {
1062     mp_result res;
1063
1064     /* Shortcut when denominators are already common */
1065     if(mp_int_compare(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) {
1066         if((res = (comb_f)(K, MP_NUMER_P(a), MP_NUMER_P(b),
1067                            MP_NUMER_P(c))) != MP_OK)
1068             return res;
1069         if((res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c))) != MP_OK)
1070             return res;
1071
1072         return s_rat_reduce(K, c);
1073     }
1074     else {
1075         mpz_t  temp[2];
1076         int    last = 0;
1077
1078         SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(a)), last);
1079         SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(b)), last);
1080
1081         if((res = mp_int_mul(K, TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK)
1082             goto CLEANUP;
1083         if((res = mp_int_mul(K, TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
1084             goto CLEANUP;
1085         if((res = (comb_f)(K, TEMP(0), TEMP(1), MP_NUMER_P(c))) != MP_OK)
1086             goto CLEANUP;
1087
1088         res = mp_int_mul(K, MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c));
1089
1090     CLEANUP:
1091         while(--last >= 0)
1092             mp_int_clear(K, TEMP(last));
1093
1094         if(res == MP_OK)
1095             return s_rat_reduce(K, c);
1096         else
1097             return res;
1098     }
1099 }
1100
1101 /* }}} */
1102
1103 /* Here there be dragons */
```