Coverage Report

Created: 2018-12-09 11:54

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