krational.c (24575B)
1 /* 2 ** krational.c 3 ** Kernel Rationals (fixrats and bigrats) 4 ** See Copyright Notice in klisp.h 5 */ 6 7 #include <stdbool.h> 8 #include <stdint.h> 9 #include <string.h> /* for code checking '/' & '.' */ 10 #include <inttypes.h> 11 #include <ctype.h> /* for to lower */ 12 #include <math.h> 13 14 #include "krational.h" 15 #include "kinteger.h" 16 #include "kobject.h" 17 #include "kstate.h" 18 #include "kmem.h" 19 #include "kgc.h" 20 21 /* used for res & temps in operations */ 22 /* NOTE: This is to be called only with already reduced values */ 23 TValue kbigrat_new(klisp_State *K, bool sign, uint32_t num, 24 uint32_t den) 25 { 26 Bigrat *new_bigrat = klispM_new(K, Bigrat); 27 28 /* header + gc_fields */ 29 klispC_link(K, (GCObject *) new_bigrat, K_TBIGRAT, 0); 30 31 /* bigrat specific fields */ 32 /* If later changed to alloc obj: 33 GC: root bigint & put dummy value to work if garbage collections 34 happens while allocating array */ 35 new_bigrat->num.single = num; 36 new_bigrat->num.digits = &(new_bigrat->num.single); 37 new_bigrat->num.alloc = 1; 38 new_bigrat->num.used = 1; 39 new_bigrat->num.sign = sign? MP_NEG : MP_ZPOS; 40 41 new_bigrat->den.single = den; 42 new_bigrat->den.digits = &(new_bigrat->den.single); 43 new_bigrat->den.alloc = 1; 44 new_bigrat->den.used = 1; 45 new_bigrat->den.sign = MP_ZPOS; 46 47 return gc2bigrat(new_bigrat); 48 } 49 50 /* assumes src is rooted */ 51 TValue kbigrat_copy(klisp_State *K, TValue src) 52 { 53 TValue copy = kbigrat_make_simple(K); 54 /* arguments are in reverse order with respect to mp_rat_copy */ 55 UNUSED(mp_rat_init_copy(K, tv2bigrat(copy), tv2bigrat(src))); 56 return copy; 57 } 58 59 /* 60 ** read/write interface 61 */ 62 /* this works for bigrats, bigints & fixints, returns true if ok */ 63 bool krational_read(klisp_State *K, char *buf, int32_t base, TValue *out, 64 char **end) 65 { 66 TValue res = kbigrat_make_simple(K); 67 krooted_tvs_push(K, res); 68 bool ret_val = (mp_rat_read_cstring(K, tv2bigrat(res), base, 69 buf, end) == MP_OK); 70 krooted_tvs_pop(K); 71 *out = kbigrat_try_integer(K, res); 72 73 /* TODO: ideally this should be incorporated in the read code */ 74 /* detect sign after '/', and / before numbers, those are allowed 75 by imrat but not in kernel */ 76 if (ret_val) { 77 char *slash = strchr(buf, '/'); 78 if (slash != NULL && (slash == 0 || 79 (*(slash+1) == '+' || *(slash+1) == '-') || 80 (*(slash-1) == '+' || *(slash-1) == '-'))) 81 ret_val = false; 82 } 83 84 return ret_val; 85 } 86 87 /* NOTE: allow decimal for use after #e */ 88 bool krational_read_decimal(klisp_State *K, char *buf, int32_t base, TValue *out, 89 char **end, bool *out_decimalp) 90 { 91 /* NOTE: in Kernel only base ten is allowed in decimal format */ 92 klisp_assert(base == 10); 93 94 char *my_end; 95 if (!end) /* always get the last char not read */ 96 end = &my_end; 97 98 TValue res = kbigrat_make_simple(K); 99 krooted_tvs_push(K, res); 100 mp_result ret_val = mp_rat_read_ustring(K, tv2bigrat(res), base, 101 buf, end); 102 103 /* check to see if there was a decimal point, will only 104 be written to out_decimalp if no error ocurr */ 105 /* TEMP: mp_rat_read_ustring does not set *end if an error occurs. 106 * Do not let memchr() read past the end of the buffer. */ 107 bool decimalp = (ret_val == MP_OK || ret_val == MP_TRUNC) 108 ? (memchr(buf, '.', *end - buf) != NULL) 109 : false; 110 111 /* handle exponents, avoid the case where the number finishes 112 in a decimal point (this isn't allowed by kernel */ 113 if (decimalp && ret_val == MP_TRUNC && *end != buf && 114 *((*end)-1) != '.') { 115 char *ebuf = *end; 116 char el = tolower(*ebuf); 117 /* NOTE: in klisp all exponent letters map to double */ 118 if (el == 'e' || el == 's' || el == 'f' || el == 'd' || el == 'l') { 119 ebuf++; 120 TValue tv_exp_exp = kbigint_make_simple(K); 121 krooted_tvs_push(K, tv_exp_exp); 122 Bigint *exp_exp = tv2bigint(tv_exp_exp); 123 /* base should be 10 */ 124 ret_val = mp_int_read_cstring(K, exp_exp, base, ebuf, end); 125 if (ret_val == MP_OK) { 126 /* IMath doesn't have general rational exponentiation, 127 so do it manually */ 128 TValue tv_iexp = kbigint_make_simple(K); 129 krooted_tvs_push(K, tv_iexp); 130 Bigint *iexp = tv2bigint(tv_iexp); 131 UNUSED(mp_int_set_value(K, iexp, 10)); 132 bool negativep = mp_int_compare_zero(exp_exp) < 0; 133 UNUSED(mp_int_abs(K, exp_exp, exp_exp)); 134 UNUSED(mp_int_expt_full(K, iexp, exp_exp, iexp)); 135 /* iexp has 10^|exp_exp| */ 136 if (negativep) { 137 TValue tv_rexp = kbigrat_make_simple(K); 138 krooted_tvs_push(K, tv_rexp); 139 Bigrat *rexp = tv2bigrat(tv_rexp); 140 /* copy reciprocal of iexp to rexp */ 141 UNUSED(mp_rat_zero(K, rexp)); 142 UNUSED(mp_rat_add_int(K, rexp, iexp, rexp)); 143 UNUSED(mp_rat_recip(K, rexp, rexp)); 144 /* now get true number */ 145 UNUSED(mp_rat_mul(K, tv2bigrat(res), rexp, 146 tv2bigrat(res))); 147 krooted_tvs_pop(K); /* rexp */ 148 } else { /* exponent positive, can just mult int */ 149 UNUSED(mp_rat_mul_int(K, tv2bigrat(res), iexp, 150 tv2bigrat(res))); 151 } 152 krooted_tvs_pop(K); /* iexp */ 153 /* fall through, ret_val remains MP_OK */ 154 } /* else, fall through, ret_val remains != MP_OK */ 155 krooted_tvs_pop(K); /* exp_exp */ 156 } 157 } 158 159 /* TODO: ideally this should be incorporated in the read code */ 160 /* TODO: if returning MP_TRUNC adjust the end pointer returned */ 161 /* detect sign after '/', or trailing '.' or starting '/' or '.'. 162 Those are allowed by imrat but not by kernel */ 163 if (ret_val == MP_OK) { 164 char *ch = strchr(buf, '/'); 165 if (ch != NULL && (ch == 0 || 166 (*(ch+1) == '+' || *(ch+1) == '-') || 167 (*(ch-1) == '+' || *(ch-1) == '-'))) { 168 ret_val = MP_TRUNC; 169 } else { 170 ch = strchr(buf, '.'); 171 if (ch != NULL && (ch == 0 || 172 (*(ch+1) == '+' || *(ch+1) == '-') || 173 (*(ch-1) == '+' || *(ch-1) == '-') || 174 ch == *end - 1)) { 175 ret_val = MP_TRUNC; 176 } 177 } 178 } 179 180 if (ret_val == MP_OK && out_decimalp != NULL) 181 *out_decimalp = decimalp; 182 183 krooted_tvs_pop(K); 184 if (ret_val == MP_OK) { 185 *out = kbigrat_try_integer(K, res); 186 } 187 188 return ret_val == MP_OK; 189 } 190 191 /* this is used by write to estimate the number of chars necessary to 192 print the number */ 193 int32_t kbigrat_print_size(TValue tv_bigrat, int32_t base) 194 { 195 klisp_assert(ttisbigrat(tv_bigrat)); 196 return mp_rat_string_len(tv2bigrat(tv_bigrat), base); 197 } 198 199 /* this is used by write */ 200 void kbigrat_print_string(klisp_State *K, TValue tv_bigrat, int32_t base, 201 char *buf, int32_t limit) 202 { 203 klisp_assert(ttisbigrat(tv_bigrat)); 204 mp_result res = mp_rat_to_string(K, tv2bigrat(tv_bigrat), base, buf, 205 limit); 206 /* only possible error is truncation */ 207 klisp_assert(res == MP_OK); 208 } 209 210 /* Interface for kgnumbers */ 211 212 /* The compare predicates take a klisp_State because in general 213 may need to do multiplications */ 214 bool kbigrat_eqp(klisp_State *K, TValue tv_bigrat1, TValue tv_bigrat2) 215 { 216 return (mp_rat_compare(K, tv2bigrat(tv_bigrat1), 217 tv2bigrat(tv_bigrat2)) == 0); 218 } 219 220 bool kbigrat_ltp(klisp_State *K, TValue tv_bigrat1, TValue tv_bigrat2) 221 { 222 return (mp_rat_compare(K, tv2bigrat(tv_bigrat1), 223 tv2bigrat(tv_bigrat2)) < 0); 224 } 225 226 bool kbigrat_lep(klisp_State *K, TValue tv_bigrat1, TValue tv_bigrat2) 227 { 228 return (mp_rat_compare(K, tv2bigrat(tv_bigrat1), 229 tv2bigrat(tv_bigrat2)) <= 0); 230 } 231 232 bool kbigrat_gtp(klisp_State *K, TValue tv_bigrat1, TValue tv_bigrat2) 233 { 234 return (mp_rat_compare(K, tv2bigrat(tv_bigrat1), 235 tv2bigrat(tv_bigrat2)) > 0); 236 } 237 238 bool kbigrat_gep(klisp_State *K, TValue tv_bigrat1, TValue tv_bigrat2) 239 { 240 return (mp_rat_compare(K, tv2bigrat(tv_bigrat1), 241 tv2bigrat(tv_bigrat2)) >= 0); 242 } 243 244 /* 245 ** GC: All of these assume the parameters are rooted 246 */ 247 TValue kbigrat_plus(klisp_State *K, TValue n1, TValue n2) 248 { 249 TValue res = kbigrat_make_simple(K); 250 krooted_tvs_push(K, res); 251 UNUSED(mp_rat_add(K, tv2bigrat(n1), tv2bigrat(n2), tv2bigrat(res))); 252 krooted_tvs_pop(K); 253 return kbigrat_try_integer(K, res); 254 } 255 256 TValue kbigrat_times(klisp_State *K, TValue n1, TValue n2) 257 { 258 TValue res = kbigrat_make_simple(K); 259 krooted_tvs_push(K, res); 260 UNUSED(mp_rat_mul(K, tv2bigrat(n1), tv2bigrat(n2), tv2bigrat(res))); 261 krooted_tvs_pop(K); 262 return kbigrat_try_integer(K, res); 263 } 264 265 TValue kbigrat_minus(klisp_State *K, TValue n1, TValue n2) 266 { 267 TValue res = kbigrat_make_simple(K); 268 krooted_tvs_push(K, res); 269 UNUSED(mp_rat_sub(K, tv2bigrat(n1), tv2bigrat(n2), tv2bigrat(res))); 270 krooted_tvs_pop(K); 271 return kbigrat_try_integer(K, res); 272 } 273 274 /* NOTE: n2 can't be zero, that case should be checked before calling this */ 275 TValue kbigrat_div_mod(klisp_State *K, TValue n1, TValue n2, TValue *res_r) 276 { 277 /* NOTE: quotient is always an integer, remainder may be any rational */ 278 TValue tv_q = kbigint_make_simple(K); 279 krooted_tvs_push(K, tv_q); 280 TValue tv_r = kbigint_make_simple(K); 281 krooted_tvs_push(K, tv_r); 282 /* for temp values */ 283 TValue tv_true_rem = kbigrat_make_simple(K); 284 krooted_tvs_push(K, tv_true_rem); 285 TValue tv_div = kbigrat_make_simple(K); 286 krooted_tvs_push(K, tv_div); 287 288 Bigrat *n = tv2bigrat(n1); 289 Bigrat *d = tv2bigrat(n2); 290 291 Bigint *q = tv2bigint(tv_q); 292 Bigint *r = tv2bigint(tv_r); 293 294 Bigrat *div = tv2bigrat(tv_div); 295 Bigrat *trem = tv2bigrat(tv_true_rem); 296 297 UNUSED(mp_rat_div(K, n, d, div)); 298 299 /* Now use the integral part as the quotient and the fractional part times 300 the divisor as the remainder, but then correct the remainder so that it's 301 always positive like in kbigint_div_and_mod */ 302 303 UNUSED(mp_int_div(K, MP_NUMER_P(div), MP_DENOM_P(div), q, r)); 304 305 /* NOTE: denom is positive, so div & q & r have the same sign */ 306 307 /* first adjust the quotient if necessary, 308 the remainder will just fall into place after this */ 309 if (mp_rat_compare_zero(n) < 0) 310 UNUSED(mp_int_add_value(K, q, mp_rat_compare_zero(d) < 0? 1 : -1, q)); 311 312 UNUSED(mp_rat_sub_int(K, div, q, trem)); 313 UNUSED(mp_rat_mul(K, trem, d, trem)); 314 315 krooted_tvs_pop(K); 316 krooted_tvs_pop(K); 317 krooted_tvs_pop(K); 318 krooted_tvs_pop(K); 319 320 *res_r = kbigrat_try_integer(K, tv_true_rem); 321 return kbigrat_try_integer(K, tv_q); 322 } 323 324 TValue kbigrat_div0_mod0(klisp_State *K, TValue n1, TValue n2, TValue *res_r) 325 { 326 /* NOTE: quotient is always an integer, remainder may be any rational */ 327 TValue tv_q = kbigint_make_simple(K); 328 krooted_tvs_push(K, tv_q); 329 TValue tv_r = kbigint_make_simple(K); 330 krooted_tvs_push(K, tv_r); 331 /* for temp values */ 332 TValue tv_true_rem = kbigrat_make_simple(K); 333 krooted_tvs_push(K, tv_true_rem); 334 TValue tv_div = kbigrat_make_simple(K); 335 krooted_tvs_push(K, tv_div); 336 337 Bigrat *n = tv2bigrat(n1); 338 Bigrat *d = tv2bigrat(n2); 339 340 Bigint *q = tv2bigint(tv_q); 341 Bigint *r = tv2bigint(tv_r); 342 343 Bigrat *div = tv2bigrat(tv_div); 344 Bigrat *trem = tv2bigrat(tv_true_rem); 345 346 UNUSED(mp_rat_div(K, n, d, div)); 347 348 /* Now use the integral part as the quotient and the fractional part times 349 the divisor as the remainder, but then correct the remainder so that it's 350 in the interval [-|d/2|, |d/2|) */ 351 352 UNUSED(mp_int_div(K, MP_NUMER_P(div), MP_DENOM_P(div), q, r)); 353 /* NOTE: denom is positive, so div & q & r have the same sign */ 354 UNUSED(mp_rat_sub_int(K, div, q, trem)); 355 UNUSED(mp_rat_mul(K, trem, d, trem)); 356 357 /* NOTE: temporarily use trem as d/2 */ 358 TValue tv_d_2 = kbigrat_make_simple(K); 359 krooted_tvs_push(K, tv_d_2); 360 Bigrat *d_2 = tv2bigrat(tv_d_2); 361 TValue m2 = i2tv(2); 362 kensure_bigint(m2); 363 UNUSED(mp_rat_div_int(K, d, tv2bigint(m2), d_2)); 364 /* adjust remainder and quotient if necessary */ 365 /* first check positive side (closed part of the interval) */ 366 mp_rat_abs(K, d_2, d_2); 367 368 /* the case analysis is like in bigint (and inverse to that of fixint) */ 369 if (mp_rat_compare(K, trem, d_2) >= 0) { 370 if (mp_rat_compare_zero(d) < 0) { 371 mp_rat_add(K, trem, d, trem); 372 mp_int_sub_value(K, q, 1, q); 373 } else { 374 mp_rat_sub(K, trem, d, trem); 375 mp_int_add_value(K, q, 1, q); 376 } 377 } else { 378 /* now check negative side (open part of the interval) */ 379 mp_rat_neg(K, d_2, d_2); 380 if (mp_rat_compare(K, trem, d_2) < 0) { 381 if (mp_rat_compare_zero(d) < 0) { 382 mp_rat_sub(K, trem, d, trem); 383 mp_int_add_value(K, q, 1, q); 384 } else { 385 mp_rat_add(K, trem, d, trem); 386 mp_int_sub_value(K, q, 1, q); 387 } 388 } 389 } 390 391 krooted_tvs_pop(K); /* d/2 */ 392 393 krooted_tvs_pop(K); 394 krooted_tvs_pop(K); 395 krooted_tvs_pop(K); 396 krooted_tvs_pop(K); 397 398 *res_r = kbigrat_try_integer(K, tv_true_rem); 399 return kbigrat_try_integer(K, tv_q); 400 } 401 402 403 TValue kbigrat_divided(klisp_State *K, TValue n1, TValue n2) 404 { 405 TValue res = kbigrat_make_simple(K); 406 krooted_tvs_push(K, res); 407 UNUSED(mp_rat_div(K, tv2bigrat(n1), tv2bigrat(n2), tv2bigrat(res))); 408 krooted_tvs_pop(K); 409 return kbigrat_try_integer(K, res); 410 } 411 412 bool kbigrat_negativep(TValue tv_bigrat) 413 { 414 return (mp_rat_compare_zero(tv2bigrat(tv_bigrat)) < 0); 415 } 416 417 bool kbigrat_positivep(TValue tv_bigrat) 418 { 419 return (mp_rat_compare_zero(tv2bigrat(tv_bigrat)) > 0); 420 } 421 422 /* GC: These assume tv_bigrat is rooted */ 423 /* needs the state to create a copy if negative */ 424 TValue kbigrat_abs(klisp_State *K, TValue tv_bigrat) 425 { 426 if (kbigrat_negativep(tv_bigrat)) { 427 TValue copy = kbigrat_make_simple(K); 428 krooted_tvs_push(K, copy); 429 UNUSED(mp_rat_abs(K, tv2bigrat(tv_bigrat), tv2bigrat(copy))); 430 krooted_tvs_pop(K); 431 /* NOTE: this can never be an integer if the parameter was a bigrat */ 432 return copy; 433 } else { 434 return tv_bigrat; 435 } 436 } 437 438 TValue kbigrat_numerator(klisp_State *K, TValue tv_bigrat) 439 { 440 int32_t fnum = 0; 441 Bigrat *bigrat = tv2bigrat(tv_bigrat); 442 if (mp_rat_to_ints(bigrat, &fnum, NULL) == MP_OK) 443 return i2tv(fnum); 444 else { 445 TValue copy = kbigint_make_simple(K); 446 krooted_tvs_push(K, copy); 447 UNUSED(mp_rat_numer(K, bigrat, tv2bigint(copy))); 448 krooted_tvs_pop(K); 449 /* NOTE: may still be a fixint because mp_rat_to_ints fails if 450 either numer or denom isn't a fixint */ 451 return kbigint_try_fixint(K, copy); 452 } 453 } 454 455 TValue kbigrat_denominator(klisp_State *K, TValue tv_bigrat) 456 { 457 int32_t fden = 0; 458 Bigrat *bigrat = tv2bigrat(tv_bigrat); 459 if (mp_rat_to_ints(bigrat, NULL, &fden) == MP_OK) 460 return i2tv(fden); 461 else { 462 TValue copy = kbigint_make_simple(K); 463 krooted_tvs_push(K, copy); 464 UNUSED(mp_rat_denom(K, bigrat, tv2bigint(copy))); 465 krooted_tvs_pop(K); 466 /* NOTE: may still be a fixint because mp_rat_to_ints fails if 467 either numer or denom isn't a fixint */ 468 return kbigint_try_fixint(K, copy); 469 } 470 } 471 472 TValue kbigrat_to_integer(klisp_State *K, TValue tv_bigrat, kround_mode mode) 473 { 474 /* do an usigned divide first */ 475 TValue tv_quot = kbigint_make_simple(K); 476 krooted_tvs_push(K, tv_quot); 477 TValue tv_rest = kbigint_make_simple(K); 478 krooted_tvs_push(K, tv_rest); 479 Bigint *quot = tv2bigint(tv_quot); 480 Bigint *rest = tv2bigint(tv_rest); 481 Bigrat *n = tv2bigrat(tv_bigrat); 482 483 UNUSED(mp_int_abs(K, MP_NUMER_P(n), quot)); 484 UNUSED(mp_int_div(K, quot, MP_DENOM_P(n), quot, rest)); 485 486 if (mp_rat_compare_zero(n) < 0) 487 UNUSED(mp_int_neg(K, quot, quot)); 488 489 switch(mode) { 490 case K_TRUNCATE: 491 /* nothing needs to be done */ 492 break; 493 case K_CEILING: 494 if (mp_rat_compare_zero(n) > 0 && mp_int_compare_zero(rest) != 0) 495 UNUSED(mp_int_add_value(K, quot, 1, quot)); 496 break; 497 case K_FLOOR: 498 if (mp_rat_compare_zero(n) < 0 && mp_int_compare_zero(rest) != 0) 499 UNUSED(mp_int_sub_value(K, quot, 1, quot)); 500 break; 501 case K_ROUND_EVEN: { 502 UNUSED(mp_int_mul_pow2(K, rest, 1, rest)); 503 int cmp = mp_int_compare(rest, MP_DENOM_P(n)); 504 if (cmp > 0 || (cmp == 0 && mp_int_is_odd(quot))) { 505 UNUSED(mp_int_add_value(K, quot, mp_rat_compare_zero(n) < 0? 506 -1 : 1, quot)); 507 } 508 break; 509 } 510 } 511 512 krooted_tvs_pop(K); 513 krooted_tvs_pop(K); 514 return kbigint_try_fixint(K, tv_quot); 515 } 516 517 /* 518 ** SOURCE NOTE: this implementation is from the Haskell 98 report 519 */ 520 /* 521 approxRational x eps = simplest (x-eps) (x+eps) 522 where simplest x y | y < x = simplest y x 523 | x == y = xr 524 | x > 0 = simplest' n d n' d' 525 | y < 0 = - simplest' (-n') d' (-n) d 526 | otherwise = 0 :% 1 527 where xr@(n:%d) = toRational x 528 (n':%d') = toRational y 529 530 simplest' n d n' d' -- assumes 0 < n%d < n'%d' 531 | r == 0 = q :% 1 532 | q /= q' = (q+1) :% 1 533 | otherwise = (q*n''+d'') :% n'' 534 where (q,r) = quotRem n d 535 (q',r') = quotRem n' d' 536 (n'':%d'') = simplest' d' r' d r 537 538 */ 539 540 /* 541 ** NOTE: this reads almost like a Haskell commercial. 542 ** The c code is an order of magnitude longer. Some of this has to do 543 ** with the representation of values, some because this is iterative, 544 ** some because of GC rooting, some because of lack of powerful pattern 545 ** matching, and so on, and so on 546 */ 547 548 /* Assumes 0 < n1/d1 < n2/d2 */ 549 /* GC: Assumes n1, d1, n2, d2, and res are fresh (can be mutated) and rooted */ 550 static void simplest(klisp_State *K, TValue tv_n1, TValue tv_d1, 551 TValue tv_n2, TValue tv_d2, TValue tv_res) 552 { 553 Bigint *n1 = tv2bigint(tv_n1); 554 Bigint *d1 = tv2bigint(tv_d1); 555 Bigint *n2 = tv2bigint(tv_n2); 556 Bigint *d2 = tv2bigint(tv_d2); 557 558 Bigrat *res = tv2bigrat(tv_res); 559 560 /* create vars q1, r1, q2 & r2 */ 561 TValue tv_q1 = kbigint_make_simple(K); 562 krooted_tvs_push(K, tv_q1); 563 Bigint *q1 = tv2bigint(tv_q1); 564 565 TValue tv_r1 = kbigint_make_simple(K); 566 krooted_tvs_push(K, tv_r1); 567 Bigint *r1 = tv2bigint(tv_r1); 568 569 TValue tv_q2 = kbigint_make_simple(K); 570 krooted_tvs_push(K, tv_q2); 571 Bigint *q2 = tv2bigint(tv_q2); 572 573 TValue tv_r2 = kbigint_make_simple(K); 574 krooted_tvs_push(K, tv_r2); 575 Bigint *r2 = tv2bigint(tv_r2); 576 577 while(true) { 578 UNUSED(mp_int_div(K, n1, d1, q1, r1)); 579 UNUSED(mp_int_div(K, n2, d2, q2, r2)); 580 581 if (mp_int_compare_zero(r1) == 0) { 582 /* res = q1 / 1 */ 583 mp_rat_zero(K, res); 584 mp_rat_add_int(K, res, q1, res); 585 break; 586 } else if (mp_int_compare(q1, q2) != 0) { 587 /* res = (q1 + 1) / 1 */ 588 mp_rat_zero(K, res); 589 mp_int_add_value(K, q1, 1, q1); 590 mp_rat_add_int(K, res, q1, res); 591 break; 592 } else { 593 /* simulate a recursive call */ 594 TValue saved_q1 = kbigint_make_simple(K); 595 krooted_tvs_push(K, saved_q1); 596 UNUSED(mp_int_copy(K, q1, tv2bigint(saved_q1))); 597 ks_spush(K, saved_q1); 598 krooted_tvs_pop(K); 599 600 UNUSED(mp_int_copy(K, d2, n1)); 601 UNUSED(mp_int_copy(K, d1, n2)); 602 UNUSED(mp_int_copy(K, r2, d1)); 603 UNUSED(mp_int_copy(K, r1, d2)); 604 } /* fall through */ 605 } 606 607 /* now, if there were "recursive" calls, complete them */ 608 while(!ks_sisempty(K)) { 609 TValue saved_q1 = ks_sget(K); 610 TValue tv_tmp = kbigrat_make_simple(K); 611 krooted_tvs_push(K, tv_tmp); 612 Bigrat *tmp = tv2bigrat(tv_tmp); 613 614 UNUSED(mp_rat_copy(K, res, tmp)); 615 /* res = (saved_q * n(res)) + d(res)) / n(res) */ 616 UNUSED(mp_rat_zero(K, res)); 617 UNUSED(mp_rat_add_int(K, res, tv2bigint(saved_q1), res)); 618 UNUSED(mp_rat_mul_int(K, res, MP_NUMER_P(tmp), res)); 619 UNUSED(mp_rat_add_int(K, res, MP_DENOM_P(tmp), res)); 620 UNUSED(mp_rat_div_int(K, res, MP_NUMER_P(tmp), res)); 621 krooted_tvs_pop(K); /* tmp */ 622 ks_sdpop(K); /* saved_q1 */ 623 } 624 625 krooted_tvs_pop(K); /* q1, r1, q2, r2 */ 626 krooted_tvs_pop(K); 627 krooted_tvs_pop(K); 628 krooted_tvs_pop(K); 629 630 return; 631 } 632 633 TValue kbigrat_simplest_rational(klisp_State *K, TValue tv_n1, TValue tv_n2) 634 { 635 TValue tv_res = kbigrat_make_simple(K); 636 krooted_tvs_push(K, tv_res); 637 Bigrat *res = tv2bigrat(tv_res); 638 Bigrat *n1 = tv2bigrat(tv_n1); 639 Bigrat *n2 = tv2bigrat(tv_n2); 640 641 int32_t cmp = mp_rat_compare(K, n1, n2); 642 if (cmp > 0) { /* n1 > n2, swap */ 643 TValue temp = tv_n1; 644 tv_n1 = tv_n2; 645 tv_n2 = temp; 646 n1 = tv2bigrat(tv_n1); 647 n2 = tv2bigrat(tv_n2); 648 /* fall through */ 649 } else if (cmp == 0) { /* n1 == n2 */ 650 krooted_tvs_pop(K); 651 return tv_n1; 652 } /* else fall through */ 653 654 /* we now know that n1 < n2 */ 655 if (mp_rat_compare_zero(n1) > 0) { 656 /* 0 > n1 > n2 */ 657 TValue num1 = kbigint_make_simple(K); 658 krooted_tvs_push(K, num1); 659 UNUSED(mp_rat_numer(K, n1, tv2bigint(num1))); 660 661 TValue den1 = kbigint_make_simple(K); 662 krooted_tvs_push(K, den1); 663 UNUSED(mp_rat_denom(K, n1, tv2bigint(den1))); 664 665 TValue num2 = kbigint_make_simple(K); 666 krooted_tvs_push(K, num2); 667 UNUSED(mp_rat_numer(K, n2, tv2bigint(num2))); 668 669 TValue den2 = kbigint_make_simple(K); 670 krooted_tvs_push(K, den2); 671 UNUSED(mp_rat_denom(K, n2, tv2bigint(den2))); 672 673 simplest(K, num1, den1, num2, den2, tv_res); 674 675 krooted_tvs_pop(K); /* num1, num2, den1, den2 */ 676 krooted_tvs_pop(K); 677 krooted_tvs_pop(K); 678 krooted_tvs_pop(K); 679 680 krooted_tvs_pop(K); /* tv_res */ 681 return kbigrat_try_integer(K, tv_res); 682 } else if (mp_rat_compare_zero(n2) < 0) { 683 /* n1 < n2 < 0 */ 684 685 /* do -(simplest -n2/d2 -n1/d1) */ 686 687 TValue num1 = kbigint_make_simple(K); 688 krooted_tvs_push(K, num1); 689 UNUSED(mp_int_neg(K, MP_NUMER_P(n2), tv2bigint(num1))); 690 691 TValue den1 = kbigint_make_simple(K); 692 krooted_tvs_push(K, den1); 693 UNUSED(mp_rat_denom(K, n2, tv2bigint(den1))); 694 695 TValue num2 = kbigint_make_simple(K); 696 krooted_tvs_push(K, num2); 697 UNUSED(mp_int_neg(K, MP_NUMER_P(n1), tv2bigint(num2))); 698 699 TValue den2 = kbigint_make_simple(K); 700 krooted_tvs_push(K, den2); 701 UNUSED(mp_rat_denom(K, n1, tv2bigint(den2))); 702 703 simplest(K, num1, den1, num2, den2, tv_res); 704 mp_rat_neg(K, res, res); 705 706 krooted_tvs_pop(K); /* num1, num2, den1, den2 */ 707 krooted_tvs_pop(K); 708 krooted_tvs_pop(K); 709 krooted_tvs_pop(K); 710 711 krooted_tvs_pop(K); /* tv_res */ 712 return kbigrat_try_integer(K, tv_res); 713 } else { 714 /* n1 < 0 < n2 */ 715 krooted_tvs_pop(K); 716 return i2tv(0); 717 } 718 } 719 720 TValue kbigrat_rationalize(klisp_State *K, TValue tv_n1, TValue tv_n2) 721 { 722 /* delegate all work to simplest_rational */ 723 TValue tv_min = kbigrat_make_simple(K); 724 krooted_tvs_push(K, tv_min); 725 TValue tv_max = kbigrat_make_simple(K); 726 krooted_tvs_push(K, tv_max); 727 728 Bigrat *n1 = tv2bigrat(tv_n1); 729 Bigrat *n2 = tv2bigrat(tv_n2); 730 /* it doesn't matter if these are reversed */ 731 Bigrat *min = tv2bigrat(tv_min); 732 Bigrat *max = tv2bigrat(tv_max); 733 734 UNUSED(mp_rat_sub(K, n1, n2, min)); 735 UNUSED(mp_rat_add(K, n1, n2, max)); 736 737 TValue res = kbigrat_simplest_rational(K, tv_min, tv_max); 738 739 krooted_tvs_pop(K); 740 krooted_tvs_pop(K); 741 742 return res; 743 }