Coverage Report

Created: 2017-04-29 12:21

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/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
61.1k
#define TEMP(K) (temp + (K))
34
15.2k
#define SETUP(E, C) \
35
15.2k
do
{if(15.2k
(res = (E)) != MP_OK15.2k
)
goto CLEANUP0
;
++(C);}15.2k
while(
015.2k
)
36
37
/* Argument checking:
38
   Use CHECK() where a return value is required; NRCHECK() elsewhere */
39
#define CHECK(TEST)   assert(TEST)
40
33.4k
#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
33.4k
{
53
33.4k
  mp_result res;
54
33.4k
55
33.4k
  if (
(res = mp_int_init(33.4k
MP_NUMER_P33.4k
(r))) != MP_OK)
56
0
    return res;
57
33.4k
  
if (33.4k
(res = mp_int_init(33.4k
MP_DENOM_P33.4k
(r))) != MP_OK)
{0
58
0
    mp_int_clear(MP_NUMER_P(r));
59
0
    return res;
60
0
  }
61
33.4k
62
33.4k
  
return mp_int_set_value(33.4k
MP_DENOM_P33.4k
(r), 1);
63
33.4k
}
64
65
mp_rat mp_rat_alloc(void)
66
33.4k
{
67
33.4k
  mp_rat out = malloc(sizeof(*out));
68
33.4k
69
33.4k
  if (
out != NULL33.4k
)
{33.4k
70
33.4k
    if (
mp_rat_init(out) != MP_OK33.4k
)
{0
71
0
      free(out);
72
0
      return NULL;
73
0
    }
74
33.4k
  }
75
33.4k
76
33.4k
  return out;
77
33.4k
}
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(0
MP_NUMER_P0
(r), n_prec)) != MP_OK)
88
0
    return res;
89
0
  
if (0
(res = mp_int_init_size(0
MP_DENOM_P0
(r), d_prec)) != MP_OK)
{0
90
0
    mp_int_clear(MP_NUMER_P(r));
91
0
    return res;
92
0
  }
93
0
  
94
0
  
return mp_int_set_value(0
MP_DENOM_P0
(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(0
MP_NUMER_P0
(r),
MP_NUMER_P0
(old))) != MP_OK)
102
0
    return res;
103
0
  
if (0
(res = mp_int_init_copy(0
MP_DENOM_P0
(r),
MP_DENOM_P0
(old))) != MP_OK)
104
0
    
mp_int_clear(0
MP_NUMER_P0
(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 (0
(res = mp_int_set_value(0
MP_NUMER_P0
(r), numer)) != MP_OK)
117
0
    return res;
118
0
  
if (0
(res = mp_int_set_value(0
MP_DENOM_P0
(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
19.5k
{
126
19.5k
  mp_result res;
127
19.5k
128
19.5k
  if (denom == 0)
129
0
    return MP_UNDEF;
130
19.5k
131
19.5k
  
if (19.5k
(res = mp_int_set_uvalue(19.5k
MP_NUMER_P19.5k
(r), numer)) != MP_OK)
132
0
    return res;
133
19.5k
  
if (19.5k
(res = mp_int_set_uvalue(19.5k
MP_DENOM_P19.5k
(r), denom)) != MP_OK)
134
0
    return res;
135
19.5k
136
19.5k
  return s_rat_reduce(r);
137
19.5k
}
138
139
void      mp_rat_clear(mp_rat r)
140
33.4k
{
141
33.4k
  mp_int_clear(MP_NUMER_P(r));
142
33.4k
  mp_int_clear(MP_DENOM_P(r));
143
33.4k
144
33.4k
}
145
146
void      mp_rat_free(mp_rat r)
147
33.4k
{
148
33.4k
  NRCHECK(r != NULL);
149
33.4k
  
150
33.4k
  if (r->num.digits != NULL)
151
33.4k
    mp_rat_clear(r);
152
33.4k
153
33.4k
  free(r);
154
33.4k
}
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
56.5k
{
163
56.5k
  return MP_NUMER_P(r);
164
56.5k
}
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
56.5k
{
174
56.5k
  return MP_DENOM_P(r);
175
56.5k
}
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
28.7k
{
184
28.7k
  mp_result res;
185
28.7k
186
28.7k
  if (
(res = mp_int_copy(28.7k
MP_NUMER_P28.7k
(a),
MP_NUMER_P28.7k
(c))) != MP_OK)
187
0
    return res;
188
28.7k
  
189
28.7k
  
res = mp_int_copy(28.7k
MP_DENOM_P28.7k
(a),
MP_DENOM_P28.7k
(c));
190
28.7k
  return res;
191
28.7k
}
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(0
MP_NUMER_P0
(a),
MP_NUMER_P0
(c))) != MP_OK)
205
0
    return res;
206
0
  
207
0
  
res = mp_int_abs(0
MP_DENOM_P0
(a),
MP_DENOM_P0
(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(0
MP_NUMER_P0
(a),
MP_NUMER_P0
(c))) != MP_OK)
216
0
    return res;
217
0
218
0
  
res = mp_int_copy(0
MP_DENOM_P0
(a),
MP_DENOM_P0
(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 (0
(res = mp_rat_copy(a, c)) != MP_OK0
)
230
0
    return res;
231
0
232
0
  
mp_int_swap(0
MP_NUMER_P0
(c),
MP_DENOM_P0
(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_SIGN0
(MP_NUMER_P(c)) =
MP_SIGN0
(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
16.7k
{
259
16.7k
  mp_result res;
260
16.7k
261
16.7k
  if (
(res = mp_int_mul(16.7k
MP_NUMER_P16.7k
(a),
MP_NUMER_P16.7k
(b),
MP_NUMER_P16.7k
(c))) != MP_OK)
262
0
    return res;
263
16.7k
264
16.7k
  
if (16.7k
mp_int_compare_zero(16.7k
MP_NUMER_P16.7k
(c)) != 0)
{16.6k
265
16.6k
    if (
(res = mp_int_mul(16.6k
MP_DENOM_P16.6k
(a),
MP_DENOM_P16.6k
(b),
MP_DENOM_P16.6k
(c))) != MP_OK)
266
0
      return res;
267
16.6k
  }
268
16.7k
269
16.7k
  return s_rat_reduce(c);
270
16.7k
}
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 (0
c == a || 0
c == b0
)
{0
280
0
    mpz_t tmp;
281
0
282
0
    if (
(res = mp_int_init(&tmp)) != MP_OK0
)
return res0
;
283
0
    
if (0
(res = mp_int_mul(0
MP_NUMER_P0
(a),
MP_DENOM_P0
(b), &tmp)) != MP_OK)
284
0
      goto CLEANUP;
285
0
    
if (0
(res = mp_int_mul(0
MP_DENOM_P0
(a),
MP_NUMER_P0
(b),
MP_DENOM_P0
(c))) != MP_OK)
286
0
      goto CLEANUP;
287
0
    
res = mp_int_copy(&tmp, 0
MP_NUMER_P0
(c));
288
0
289
0
  CLEANUP:
290
0
    mp_int_clear(&tmp);
291
0
  }
292
0
  else {
293
0
    if (
(res = mp_int_mul(0
MP_NUMER_P0
(a),
MP_DENOM_P0
(b),
MP_NUMER_P0
(c))) != MP_OK)
294
0
      return res;
295
0
    
if (0
(res = mp_int_mul(0
MP_DENOM_P0
(a),
MP_NUMER_P0
(b),
MP_DENOM_P0
(c))) != MP_OK)
296
0
      return res;
297
0
  }
298
0
299
0
  
if (0
res != MP_OK0
)
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 (0
(res = mp_int_mul(&tmp, 0
MP_DENOM_P0
(a), &tmp)) != MP_OK)
314
0
    goto CLEANUP;
315
0
316
0
  
if (0
(res = mp_rat_copy(a, c)) != MP_OK0
)
317
0
    goto CLEANUP;
318
0
319
0
  
if (0
(res = mp_int_add(0
MP_NUMER_P0
(c), &tmp,
MP_NUMER_P0
(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 (0
(res = mp_int_mul(&tmp, 0
MP_DENOM_P0
(a), &tmp)) != MP_OK)
338
0
    goto CLEANUP;
339
0
340
0
  
if (0
(res = mp_rat_copy(a, c)) != MP_OK0
)
341
0
    goto CLEANUP;
342
0
343
0
  
if (0
(res = mp_int_sub(0
MP_NUMER_P0
(c), &tmp,
MP_NUMER_P0
(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 (0
(res = mp_int_mul(0
MP_NUMER_P0
(c), b,
MP_NUMER_P0
(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 (0
(res = mp_rat_copy(a, c)) != MP_OK0
)
374
0
    return res;
375
0
376
0
  
if (0
(res = mp_int_mul(0
MP_DENOM_P0
(c), b,
MP_DENOM_P0
(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(0
b == 10
)
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 (0
(res = mp_int_expt(0
MP_NUMER_P0
(a), b,
MP_NUMER_P0
(c))) != MP_OK)
395
0
    return res;
396
0
397
0
  
return mp_int_expt(0
MP_DENOM_P0
(a), b,
MP_DENOM_P0
(c));
398
0
}
399
400
int       mp_rat_compare(mp_rat a, mp_rat b)
401
13.7k
{
402
13.7k
  /* Quick check for opposite signs.  Works because the sign of the numerator
403
13.7k
     is always definitive. */
404
13.7k
  if (
MP_SIGN13.7k
(MP_NUMER_P(a)) != 13.7k
MP_SIGN13.7k
(MP_NUMER_P(b)))
{0
405
0
    if (
MP_SIGN0
(MP_NUMER_P(a)) == MP_ZPOS0
)
406
0
      return 1;
407
0
    else
408
0
      return -1;
409
0
  }
410
13.7k
  else {
411
13.7k
    /* Compare absolute magnitudes; if both are positive, the answer stands,
412
13.7k
       otherwise it needs to be reflected about zero. */
413
13.7k
    int cmp = mp_rat_compare_unsigned(a, b);
414
13.7k
415
13.7k
    if (
MP_SIGN13.7k
(MP_NUMER_P(a)) == MP_ZPOS13.7k
)
416
13.7k
      return cmp;
417
13.7k
    else
418
0
      return -cmp;
419
13.7k
  }
420
13.7k
}
421
422
int       mp_rat_compare_unsigned(mp_rat a, mp_rat b)
423
13.7k
{
424
13.7k
  /* If the denominators are equal, we can quickly compare numerators without
425
13.7k
     multiplying.  Otherwise, we actually have to do some work. */
426
13.7k
  if (
mp_int_compare_unsigned(13.7k
MP_DENOM_P13.7k
(a),
MP_DENOM_P13.7k
(b)) == 0)
427
6.09k
    
return mp_int_compare_unsigned(6.09k
MP_NUMER_P6.09k
(a),
MP_NUMER_P6.09k
(b));
428
13.7k
429
7.64k
  else {
430
7.64k
    mpz_t  temp[2];
431
7.64k
    mp_result res;
432
7.64k
    int  cmp = INT_MAX, last = 0;
433
7.64k
434
7.64k
    /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
435
7.64k
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
436
7.64k
    
SETUP7.64k
(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);7.64k
437
7.64k
438
7.64k
    
if (7.64k
(res = mp_int_mul(7.64k
TEMP7.64k
(0),
MP_DENOM_P7.64k
(b),
TEMP7.64k
(0))) != MP_OK ||
439
7.64k
  
(res = mp_int_mul(7.64k
TEMP7.64k
(1),
MP_DENOM_P7.64k
(a),
TEMP7.64k
(1))) != MP_OK)
440
0
      goto CLEANUP;
441
7.64k
    
442
7.64k
    
cmp = mp_int_compare_unsigned(7.64k
TEMP7.64k
(0),
TEMP7.64k
(1));
443
7.64k
    
444
7.64k
  CLEANUP:
445
22.9k
    while (--last >= 0)
446
15.2k
      
mp_int_clear(15.2k
TEMP15.2k
(last));
447
7.64k
448
7.64k
    return cmp;
449
7.64k
  }
450
13.7k
}
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 (0
(res = mp_rat_set_value(&tmp, n, d)) != MP_OK0
)
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(0
MP_NUMER_P0
(r), num)) != MP_OK)
485
0
    return res;
486
0
487
0
  
res = mp_int_to_int(0
MP_DENOM_P0
(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(0
MP_NUMER_P0
(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 (0
mp_int_compare_zero(0
MP_NUMER_P0
(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
  
SETUP0
(mp_int_init(TEMP(last)), last);0
531
0
  
SETUP0
(mp_int_init(TEMP(last)), last);0
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(0
TEMP0
(0),
TEMP0
(0));
536
0
  if (
(res = mp_int_div(0
TEMP0
(0),
MP_DENOM_P0
(r),
TEMP0
(0),
TEMP0
(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; 0
lead_0 < prec && 0
mp_int_compare(0
TEMP0
(1),
MP_DENOM_P0
(r)) < 0;
544
0
      
++lead_00
)
{0
545
0
    if (
(res = mp_int_mul_value(0
TEMP0
(1), radix,
TEMP0
(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 (0
prec > lead_00
)
{0
552
0
    if (
(res = mp_int_expt_value(radix, prec - lead_0, 0
TEMP0
(2))) != MP_OK)
553
0
      goto CLEANUP;
554
0
    
if (0
(res = mp_int_mul(0
TEMP0
(1),
TEMP0
(2),
TEMP0
(1))) != MP_OK)
555
0
      goto CLEANUP;
556
0
  }
557
0
  
if (0
(res = mp_int_div(0
TEMP0
(1),
MP_DENOM_P0
(r),
TEMP0
(1),
TEMP0
(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(0
TEMP0
(2)) != 0)
{0
571
0
      if (prec == 0)
572
0
  
res = mp_int_add_value(0
TEMP0
(0), 1,
TEMP0
(0));
573
0
      else
574
0
  
res = mp_int_add_value(0
TEMP0
(1), 1,
TEMP0
(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(0
TEMP0
(2), 1,
TEMP0
(2))) != MP_OK)
581
0
      goto CLEANUP;
582
0
583
0
    
cmp = mp_int_compare(0
TEMP0
(2),
MP_DENOM_P0
(r));
584
0
585
0
    if (round == MP_ROUND_HALF_UP)
586
0
      cmp += 1;
587
0
588
0
    if (
cmp > 00
)
{0
589
0
      if (prec == 0)
590
0
  
res = mp_int_add_value(0
TEMP0
(0), 1,
TEMP0
(0));
591
0
      else
592
0
  
res = mp_int_add_value(0
TEMP0
(1), 1,
TEMP0
(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 (0
MP_SIGN0
(MP_NUMER_P(r)) == MP_NEG &&0
607
0
      
(mp_int_compare_zero(0
TEMP0
(0)) != 0 ||
608
0
       
mp_int_compare_zero(0
TEMP0
(1)) != 0))
{0
609
0
    *start++ = '-';
610
0
    left -= 1;
611
0
  }
612
0
613
0
  if (
(res = mp_int_to_string(0
TEMP0
(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 + 10
)
{0
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(0
TEMP0
(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(0
MP_NUMER_P0
(r)) != 0)
651
0
    
d_len = mp_int_string_len(0
MP_DENOM_P0
(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(0
MP_NUMER_P0
(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 (0
isspace((unsigned char) *endp)0
)
693
0
    ++endp;
694
0
  
695
0
  /* If there is no separator, we will stop reading at this point. */
696
0
  if (
*endp != '/'0
)
{0
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(0
MP_DENOM_P0
(r), radix, endp, end)) != MP_OK)
705
0
    return res;
706
0
  
707
0
  /* Make sure the value is well-defined */
708
0
  
if (0
mp_int_compare_zero(0
MP_DENOM_P0
(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_OK0
)
{0
730
0
    if (
res == MP_TRUNC0
)
{0
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 (0
end != NULL0
)
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 (0
(res = mp_int_read_cstring(0
MP_NUMER_P0
(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(0
MP_DENOM_P0
(r), 1);
773
0
774
0
  if (
*endp != '.'0
)
{0
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'0
)
{0
791
0
    if (end != NULL)
792
0
      *end = endp;
793
0
    return MP_OK;
794
0
  }
795
0
  else 
if(0
isspace((unsigned char) *endp) || 0
*endp == '-'0
||
*endp == '+'0
)
{0
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 (0
(res = mp_int_read_cstring(&frac, radix, endp, &endp)) != MP_OK &&0
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 (0
mp_int_divisible_value(&frac, radix)0
)
820
0
      
if (0
(res = mp_int_div_value(&frac, radix, &frac, NULL)) != MP_OK0
)
821
0
  goto CLEANUP;
822
0
    
823
0
    /* Count leading zeros after the decimal point */
824
0
    
while (0
save[num_lz] == '0'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(0
MP_DENOM_P0
(r), &frac) < 0)
{0
832
0
      if (
(res = mp_int_mul_value(0
MP_DENOM_P0
(r), radix,
MP_DENOM_P0
(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 (0
num_lz > 00
)
{0
838
0
      if (
(res = mp_int_mul_value(0
MP_DENOM_P0
(r), radix,
MP_DENOM_P0
(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 (0
(res = mp_int_mul(0
MP_NUMER_P0
(r),
MP_DENOM_P0
(r),
MP_NUMER_P0
(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(0
MP_NUMER_P0
(r), &frac,
MP_NUMER_P0
(r))) != MP_OK)
853
0
  goto CLEANUP;
854
0
855
0
      
MP_SIGN0
(MP_NUMER_P(r)) = osign;0
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
36.2k
{
880
36.2k
  mpz_t gcd;
881
36.2k
  mp_result res = MP_OK;
882
36.2k
883
36.2k
  if (
mp_int_compare_zero(36.2k
MP_NUMER_P36.2k
(r)) == 0)
{128
884
128
    mp_int_set_value(MP_DENOM_P(r), 1);
885
128
    return MP_OK;
886
128
  }
887
36.2k
888
36.2k
  /* If the greatest common divisor of the numerator and denominator is greater
889
36.2k
     than 1, divide it out. */
890
36.1k
  
if (36.1k
(res = mp_int_init(&gcd)) != MP_OK36.1k
)
891
0
    return res;
892
36.1k
893
36.1k
  
if (36.1k
(res = mp_int_gcd(36.1k
MP_NUMER_P36.1k
(r),
MP_DENOM_P36.1k
(r), &gcd)) != MP_OK)
894
0
    goto CLEANUP;
895
36.1k
896
36.1k
  
if (36.1k
mp_int_compare_value(&gcd, 1) != 036.1k
)
{8.71k
897
8.71k
    if (
(res = mp_int_div(8.71k
MP_NUMER_P8.71k
(r), &gcd,
MP_NUMER_P8.71k
(r), NULL)) != MP_OK)
898
0
      goto CLEANUP;
899
8.71k
    
if (8.71k
(res = mp_int_div(8.71k
MP_DENOM_P8.71k
(r), &gcd,
MP_DENOM_P8.71k
(r), NULL)) != MP_OK)
900
0
      goto CLEANUP;
901
8.71k
  }
902
36.1k
903
36.1k
  /* Fix up the signs of numerator and denominator */
904
36.1k
  
if (36.1k
MP_SIGN36.1k
(MP_NUMER_P(r)) == 36.1k
MP_SIGN36.1k
(MP_DENOM_P(r)))
905
36.1k
    
MP_SIGN36.1k
(MP_NUMER_P(r)) = 36.1k
MP_SIGN36.1k
(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
36.1k
911
36.1k
 CLEANUP:
912
36.1k
  mp_int_clear(&gcd);
913
36.1k
914
36.1k
  return res;
915
36.1k
}
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(0
MP_DENOM_P0
(a),
MP_DENOM_P0
(b)) == 0)
{0
924
0
    if (
(res = (comb_f)(0
MP_NUMER_P0
(a),
MP_NUMER_P0
(b),
MP_NUMER_P0
(c))) != MP_OK)
925
0
      return res;
926
0
    
if (0
(res = mp_int_copy(0
MP_DENOM_P0
(a),
MP_DENOM_P0
(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
    
SETUP0
(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);0
937
0
    
938
0
    
if (0
(res = mp_int_mul(0
TEMP0
(0),
MP_DENOM_P0
(b),
TEMP0
(0))) != MP_OK)
939
0
      goto CLEANUP;
940
0
    
if (0
(res = mp_int_mul(0
TEMP0
(1),
MP_DENOM_P0
(a),
TEMP0
(1))) != MP_OK)
941
0
      goto CLEANUP;
942
0
    
if (0
(res = (comb_f)(0
TEMP0
(0),
TEMP0
(1),
MP_NUMER_P0
(c))) != MP_OK)
943
0
      goto CLEANUP;
944
0
945
0
    
res = mp_int_mul(0
MP_DENOM_P0
(a),
MP_DENOM_P0
(b),
MP_DENOM_P0
(c));
946
0
947
0
  CLEANUP:
948
0
    while (--last >= 0) 
949
0
      
mp_int_clear(0
TEMP0
(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 */