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