Coverage Report

Created: 2017-06-23 12:40

/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
196k
#define TEMP(K) (temp + (K))
34
49.0k
#define SETUP(E, C) \
35
49.0k
do
{if(49.0k
(res = (E)) != MP_OK49.0k
)
goto CLEANUP0
;
++(C);}49.0k
while(
049.0k
)
36
37
/* Argument checking:
38
   Use CHECK() where a return value is required; NRCHECK() elsewhere */
39
#define CHECK(TEST)   assert(TEST)
40
79.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
79.4k
{
53
79.4k
  mp_result res;
54
79.4k
55
79.4k
  if (
(res = mp_int_init(79.4k
MP_NUMER_P79.4k
(r))) != MP_OK)
56
0
    return res;
57
79.4k
  
if (79.4k
(res = mp_int_init(79.4k
MP_DENOM_P79.4k
(r))) != MP_OK)
{0
58
0
    mp_int_clear(MP_NUMER_P(r));
59
0
    return res;
60
0
  }
61
79.4k
62
79.4k
  
return mp_int_set_value(79.4k
MP_DENOM_P79.4k
(r), 1);
63
79.4k
}
64
65
mp_rat mp_rat_alloc(void)
66
79.4k
{
67
79.4k
  mp_rat out = malloc(sizeof(*out));
68
79.4k
69
79.4k
  if (
out != NULL79.4k
)
{79.4k
70
79.4k
    if (
mp_rat_init(out) != MP_OK79.4k
)
{0
71
0
      free(out);
72
0
      return NULL;
73
0
    }
74
79.4k
  }
75
79.4k
76
79.4k
  return out;
77
79.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
45.3k
{
126
45.3k
  mp_result res;
127
45.3k
128
45.3k
  if (denom == 0)
129
0
    return MP_UNDEF;
130
45.3k
131
45.3k
  
if (45.3k
(res = mp_int_set_uvalue(45.3k
MP_NUMER_P45.3k
(r), numer)) != MP_OK)
132
0
    return res;
133
45.3k
  
if (45.3k
(res = mp_int_set_uvalue(45.3k
MP_DENOM_P45.3k
(r), denom)) != MP_OK)
134
0
    return res;
135
45.3k
136
45.3k
  return s_rat_reduce(r);
137
45.3k
}
138
139
void      mp_rat_clear(mp_rat r)
140
79.4k
{
141
79.4k
  mp_int_clear(MP_NUMER_P(r));
142
79.4k
  mp_int_clear(MP_DENOM_P(r));
143
79.4k
144
79.4k
}
145
146
void      mp_rat_free(mp_rat r)
147
79.4k
{
148
79.4k
  NRCHECK(r != NULL);
149
79.4k
  
150
79.4k
  if (r->num.digits != NULL)
151
79.4k
    mp_rat_clear(r);
152
79.4k
153
79.4k
  free(r);
154
79.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
133k
{
163
133k
  return MP_NUMER_P(r);
164
133k
}
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
133k
{
174
133k
  return MP_DENOM_P(r);
175
133k
}
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
70.3k
{
184
70.3k
  mp_result res;
185
70.3k
186
70.3k
  if (
(res = mp_int_copy(70.3k
MP_NUMER_P70.3k
(a),
MP_NUMER_P70.3k
(c))) != MP_OK)
187
0
    return res;
188
70.3k
  
189
70.3k
  
res = mp_int_copy(70.3k
MP_DENOM_P70.3k
(a),
MP_DENOM_P70.3k
(c));
190
70.3k
  return res;
191
70.3k
}
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
38.7k
{
259
38.7k
  mp_result res;
260
38.7k
261
38.7k
  if (
(res = mp_int_mul(38.7k
MP_NUMER_P38.7k
(a),
MP_NUMER_P38.7k
(b),
MP_NUMER_P38.7k
(c))) != MP_OK)
262
0
    return res;
263
38.7k
264
38.7k
  
if (38.7k
mp_int_compare_zero(38.7k
MP_NUMER_P38.7k
(c)) != 0)
{38.6k
265
38.6k
    if (
(res = mp_int_mul(38.6k
MP_DENOM_P38.6k
(a),
MP_DENOM_P38.6k
(b),
MP_DENOM_P38.6k
(c))) != MP_OK)
266
0
      return res;
267
38.6k
  }
268
38.7k
269
38.7k
  return s_rat_reduce(c);
270
38.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
35.3k
{
402
35.3k
  /* Quick check for opposite signs.  Works because the sign of the numerator
403
35.3k
     is always definitive. */
404
35.3k
  if (
MP_SIGN35.3k
(MP_NUMER_P(a)) != 35.3k
MP_SIGN35.3k
(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
35.3k
  else {
411
35.3k
    /* Compare absolute magnitudes; if both are positive, the answer stands,
412
35.3k
       otherwise it needs to be reflected about zero. */
413
35.3k
    int cmp = mp_rat_compare_unsigned(a, b);
414
35.3k
415
35.3k
    if (
MP_SIGN35.3k
(MP_NUMER_P(a)) == MP_ZPOS35.3k
)
416
35.3k
      return cmp;
417
35.3k
    else
418
0
      return -cmp;
419
35.3k
  }
420
35.3k
}
421
422
int       mp_rat_compare_unsigned(mp_rat a, mp_rat b)
423
35.3k
{
424
35.3k
  /* If the denominators are equal, we can quickly compare numerators without
425
35.3k
     multiplying.  Otherwise, we actually have to do some work. */
426
35.3k
  if (
mp_int_compare_unsigned(35.3k
MP_DENOM_P35.3k
(a),
MP_DENOM_P35.3k
(b)) == 0)
427
10.8k
    
return mp_int_compare_unsigned(10.8k
MP_NUMER_P10.8k
(a),
MP_NUMER_P10.8k
(b));
428
35.3k
429
24.5k
  else {
430
24.5k
    mpz_t  temp[2];
431
24.5k
    mp_result res;
432
24.5k
    int  cmp = INT_MAX, last = 0;
433
24.5k
434
24.5k
    /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
435
24.5k
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
436
24.5k
    
SETUP24.5k
(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);24.5k
437
24.5k
438
24.5k
    
if (24.5k
(res = mp_int_mul(24.5k
TEMP24.5k
(0),
MP_DENOM_P24.5k
(b),
TEMP24.5k
(0))) != MP_OK ||
439
24.5k
  
(res = mp_int_mul(24.5k
TEMP24.5k
(1),
MP_DENOM_P24.5k
(a),
TEMP24.5k
(1))) != MP_OK)
440
0
      goto CLEANUP;
441
24.5k
    
442
24.5k
    
cmp = mp_int_compare_unsigned(24.5k
TEMP24.5k
(0),
TEMP24.5k
(1));
443
24.5k
    
444
24.5k
  CLEANUP:
445
73.5k
    while (--last >= 0)
446
49.0k
      
mp_int_clear(49.0k
TEMP49.0k
(last));
447
24.5k
448
24.5k
    return cmp;
449
24.5k
  }
450
35.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 (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
84.1k
{
880
84.1k
  mpz_t gcd;
881
84.1k
  mp_result res = MP_OK;
882
84.1k
883
84.1k
  if (
mp_int_compare_zero(84.1k
MP_NUMER_P84.1k
(r)) == 0)
{161
884
161
    mp_int_set_value(MP_DENOM_P(r), 1);
885
161
    return MP_OK;
886
161
  }
887
84.1k
888
84.1k
  /* If the greatest common divisor of the numerator and denominator is greater
889
84.1k
     than 1, divide it out. */
890
83.9k
  
if (83.9k
(res = mp_int_init(&gcd)) != MP_OK83.9k
)
891
0
    return res;
892
83.9k
893
83.9k
  
if (83.9k
(res = mp_int_gcd(83.9k
MP_NUMER_P83.9k
(r),
MP_DENOM_P83.9k
(r), &gcd)) != MP_OK)
894
0
    goto CLEANUP;
895
83.9k
896
83.9k
  
if (83.9k
mp_int_compare_value(&gcd, 1) != 083.9k
)
{22.5k
897
22.5k
    if (
(res = mp_int_div(22.5k
MP_NUMER_P22.5k
(r), &gcd,
MP_NUMER_P22.5k
(r), NULL)) != MP_OK)
898
0
      goto CLEANUP;
899
22.5k
    
if (22.5k
(res = mp_int_div(22.5k
MP_DENOM_P22.5k
(r), &gcd,
MP_DENOM_P22.5k
(r), NULL)) != MP_OK)
900
0
      goto CLEANUP;
901
22.5k
  }
902
83.9k
903
83.9k
  /* Fix up the signs of numerator and denominator */
904
83.9k
  
if (83.9k
MP_SIGN83.9k
(MP_NUMER_P(r)) == 83.9k
MP_SIGN83.9k
(MP_DENOM_P(r)))
905
83.9k
    
MP_SIGN83.9k
(MP_NUMER_P(r)) = 83.9k
MP_SIGN83.9k
(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
83.9k
911
83.9k
 CLEANUP:
912
83.9k
  mp_int_clear(&gcd);
913
83.9k
914
83.9k
  return res;
915
83.9k
}
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 */