kreal.c (20981B)
1 /* 2 ** kreal.c 3 ** Kernel Reals (doubles) 4 ** See Copyright Notice in klisp.h 5 */ 6 7 #include <stdbool.h> 8 #include <stdint.h> 9 #include <string.h> 10 #include <inttypes.h> 11 #include <ctype.h> 12 #include <math.h> 13 #include <fenv.h> /* for setting round direction */ 14 15 #include "kreal.h" 16 #include "krational.h" 17 #include "kinteger.h" 18 #include "kobject.h" 19 #include "kstate.h" 20 #include "kmem.h" 21 #include "kgc.h" 22 #include "kpair.h" /* for list in throw error */ 23 #include "kerror.h" 24 25 double kbigint_to_double(Bigint *bigint) 26 { 27 double radix = (double) UINT32_MAX + 1.0; 28 uint32_t ndigits = bigint->used; 29 double accum = 0.0; 30 31 /* bigint is in little endian format, but we traverse in 32 big endian */ 33 do { 34 --ndigits; 35 accum = accum * radix + (double) bigint->digits[ndigits]; 36 } while (ndigits > 0); /* have to compare like this, it's unsigned */ 37 return mp_int_compare_zero(bigint) < 0? -accum : accum; 38 } 39 40 /* bigrat is rooted */ 41 double kbigrat_to_double(klisp_State *K, Bigrat *bigrat) 42 { 43 /* TODO: check rounding in extreme cases */ 44 TValue tv_rem = kbigrat_copy(K, gc2bigrat(bigrat)); 45 krooted_tvs_push(K, tv_rem); 46 Bigrat *rem = tv2bigrat(tv_rem); 47 UNUSED(mp_rat_abs(K, rem, rem)); 48 49 TValue int_radix = kbigint_make_simple(K); 50 krooted_tvs_push(K, int_radix); 51 /* cant do UINT32_MAX and then +1 because _value functions take 52 int32_t arguments */ 53 UNUSED(mp_int_set_value(K, tv2bigint(int_radix), INT32_MAX)); 54 UNUSED(mp_int_add_value(K, tv2bigint(int_radix), 1, 55 tv2bigint(int_radix))); 56 UNUSED(mp_int_add(K, tv2bigint(int_radix), tv2bigint(int_radix), 57 tv2bigint(int_radix))); 58 59 TValue int_part = kbigint_make_simple(K); 60 krooted_tvs_push(K, int_part); 61 62 double accum = 0.0; 63 double radix = (double) UINT32_MAX + 1.0; 64 uint32_t digit = 0; 65 /* inside there is a check to avoid problem with continuing fractions */ 66 while(mp_rat_compare_zero(rem) > 0) { 67 UNUSED(mp_int_div(K, MP_NUMER_P(rem), MP_DENOM_P(rem), 68 tv2bigint(int_part), NULL)); 69 70 double di = kbigint_to_double(tv2bigint(int_part)); 71 bool was_zero = di == 0.0; 72 for (uint32_t i = 0; i < digit; i++) { 73 di /= radix; 74 } 75 /* if last di became 0.0 we will exit next loop, 76 this is to avoid problem with continuing fractions */ 77 if (!was_zero && di == 0.0) 78 break; 79 80 ++digit; 81 accum += di; 82 83 UNUSED(mp_rat_sub_int(K, rem, tv2bigint(int_part), rem)); 84 UNUSED(mp_rat_mul_int(K, rem, tv2bigint(int_radix), rem)); 85 } 86 krooted_tvs_pop(K); /* int_part */ 87 krooted_tvs_pop(K); /* int_radix */ 88 krooted_tvs_pop(K); /* rem */ 89 90 return mp_rat_compare_zero(bigrat) < 0? -accum : accum; 91 } 92 93 /* TODO test strict arithmetic and throw errors on overflow and underflow? 94 if set */ 95 TValue kexact_to_inexact(klisp_State *K, TValue n) 96 { 97 bool strictp = kcurr_strict_arithp(K); 98 99 switch(ttype(n)) { 100 case K_TFIXINT: 101 /* NOTE: can't over or underflow, and can't give NaN */ 102 return d2tv((double) ivalue(n)); 103 case K_TBIGINT: { 104 Bigint *bigint = tv2bigint(n); 105 double d = kbigint_to_double(bigint); 106 if (strictp && (d == 0.0 || isinf(d) || isnan(d))) { 107 /* NOTE: bigints can't be zero */ 108 char *msg; 109 if (isnan(d)) 110 msg = "unexpected error"; 111 else if (isinf(d)) 112 msg = "overflow"; 113 else 114 msg = "undeflow"; 115 klispE_throw_simple_with_irritants(K, msg, 1, n); 116 return KUNDEF; 117 } else { 118 /* d may be inf, ktag_double will handle it */ 119 return ktag_double(d); 120 } 121 } 122 case K_TBIGRAT: { 123 Bigrat *bigrat = tv2bigrat(n); 124 double d = kbigrat_to_double(K, bigrat); 125 /* REFACTOR: this code is the same for bigints... */ 126 if (strictp && (d == 0.0 || isinf(d) || isnan(d))) { 127 /* NOTE: bigrats can't be zero */ 128 char *msg; 129 if (isnan(d)) 130 msg = "unexpected error"; 131 else if (isinf(d)) 132 msg = "overflow"; 133 else 134 msg = "undeflow"; 135 klispE_throw_simple_with_irritants(K, msg, 1, n); 136 return KUNDEF; 137 } else { 138 /* d may be inf, ktag_double will handle it */ 139 return ktag_double(d); 140 } 141 } 142 case K_TEINF: 143 return tv_equal(n, KEPINF)? KIPINF : KIMINF; 144 /* all of these are already inexact */ 145 case K_TDOUBLE: 146 case K_TIINF: 147 case K_TRWNPV: 148 case K_TUNDEFINED: 149 return n; 150 default: 151 klisp_assert(0); 152 return KUNDEF; 153 } 154 } 155 156 /* assume d is integer and doesn't fit in a fixint */ 157 TValue kdouble_to_bigint(klisp_State *K, double d) 158 { 159 bool neg = d < 0; 160 if (neg) 161 d = -d; 162 163 TValue tv_res = kbigint_make_simple(K); 164 krooted_tvs_push(K, tv_res); 165 Bigint *res = tv2bigint(tv_res); 166 mp_int_set_value(K, res, 0); 167 168 TValue tv_digit = kbigint_make_simple(K); 169 krooted_tvs_push(K, tv_digit); 170 Bigint *digit = tv2bigint(tv_digit); 171 172 /* do it 32 bits at a time */ 173 double radix = ((double) UINT32_MAX) + 1.0; 174 int power = 0; 175 while(d > 0) { 176 double dd = fmod(d, radix); 177 d = floor(d / radix); 178 /* load in two moves because set_value takes signed ints */ 179 uint32_t id = (uint32_t) dd; 180 int32_t id1 = (int32_t) (id >> 1); 181 int32_t id2 = (int32_t) (id - id1); 182 183 mp_int_set_value(K, digit, id1); 184 mp_int_add_value(K, digit, id2, digit); 185 186 mp_int_mul_pow2(K, digit, power, digit); 187 mp_int_add(K, res, digit, res); 188 189 power += 32; 190 } 191 192 if (neg) 193 mp_int_neg(K, res, res); 194 195 krooted_tvs_pop(K); /* digit */ 196 krooted_tvs_pop(K); /* res */ 197 198 return tv_res; /* can't be fixint except when coming from 199 kdouble_to_bigrat, so don't convert */ 200 } 201 202 /* TODO: should do something like rationalize with range +/- 1/2ulp) */ 203 TValue kdouble_to_bigrat(klisp_State *K, double d) 204 { 205 bool neg = d < 0; 206 if (neg) 207 d = -d; 208 209 /* find an integer, convert it and divide by 210 an adequate power of 2 */ 211 TValue tv_res = kbigrat_make_simple(K); 212 krooted_tvs_push(K, tv_res); 213 Bigrat *res = tv2bigrat(tv_res); 214 UNUSED(mp_rat_set_value(K, res, 0, 1)); 215 216 TValue tv_den = kbigint_make_simple(K); 217 krooted_tvs_push(K, tv_den); 218 Bigint *den = tv2bigint(tv_den); 219 UNUSED(mp_int_set_value(K, den, 1)); 220 221 /* XXX could be made a lot more efficiently reading ieee 222 fields directly */ 223 int ie; 224 d = frexp(d, &ie); 225 226 while(d != floor(d)) { 227 d *= 2.0; 228 --ie; 229 } 230 UNUSED(mp_int_mul_pow2(K, den, -ie, den)); 231 232 TValue tv_num = kdouble_to_bigint(K, d); 233 krooted_tvs_push(K, tv_num); 234 Bigint *num = tv2bigint(tv_num); 235 236 TValue tv_den2 = kbigrat_make_simple(K); 237 krooted_tvs_push(K, tv_den2); 238 Bigrat *den2 = tv2bigrat(tv_den2); 239 240 UNUSED(mp_rat_set_value(K, den2, 0, 1)); 241 UNUSED(mp_rat_add_int(K, den2, den, den2)); 242 UNUSED(mp_rat_set_value(K, res, 0, 1)); 243 UNUSED(mp_rat_add_int(K, res, num, res)); 244 UNUSED(mp_rat_div(K, res, den2, res)); 245 246 if (neg) 247 UNUSED(mp_rat_neg(K, res, res)); 248 249 /* now create a value corresponding to 1/2 ulp 250 for rationalize */ 251 UNUSED(mp_int_mul_pow2(K, den, 1, den)); 252 UNUSED(mp_rat_set_value(K, den2, 0, 1)); 253 UNUSED(mp_rat_add_int(K, den2, den, den2)); 254 UNUSED(mp_rat_recip(K, den2, den2)); 255 /* den2 now has 1/2 ulp */ 256 257 TValue rationalized = kbigrat_rationalize(K, tv_res, tv_den2); 258 259 krooted_tvs_pop(K); /* den2 */ 260 krooted_tvs_pop(K); /* num */ 261 krooted_tvs_pop(K); /* den */ 262 krooted_tvs_pop(K); /* res */ 263 264 return rationalized; 265 } 266 267 TValue kinexact_to_exact(klisp_State *K, TValue n) 268 { 269 switch(ttype(n)) { 270 case K_TFIXINT: 271 case K_TBIGINT: 272 case K_TBIGRAT: 273 case K_TEINF: 274 /* all of these are already exact */ 275 return n; 276 case K_TDOUBLE: { 277 double d = dvalue(n); 278 klisp_assert(!isnan(d) && !isinf(d)); 279 if (d == floor(d)) { /* integer */ 280 if (d <= (double) INT32_MAX && 281 d >= (double) INT32_MIN) { 282 return i2tv((int32_t) d); /* fixint */ 283 } else { 284 return kdouble_to_bigint(K, d); 285 } 286 } else { /* non integer rational */ 287 return kdouble_to_bigrat(K, d); 288 } 289 } 290 case K_TIINF: 291 return tv_equal(n, KIPINF)? KEPINF : KEMINF; 292 case K_TRWNPV: 293 case K_TUNDEFINED: 294 klispE_throw_simple_with_irritants(K, "no primary value", 1, n); 295 return KUNDEF; 296 default: 297 klisp_assert(0); 298 return KUNDEF; 299 } 300 } 301 302 /* 303 ** read/write interface 304 */ 305 306 /* 307 ** SOURCE NOTE: This is a more or less vanilla implementation of the algorithm 308 ** described in "How to Print Floating-Point Numbers Accurately" by 309 ** Guy L. Steele Jr. & John L. White. 310 */ 311 312 /* 313 ** TODO add awareness of read rounding (e.g. problem with 1.0e23) 314 ** TODO add exponent if too small or too big 315 */ 316 317 mp_result shift_2(klisp_State *K, Bigint *x, Bigint *n, Bigint *r) 318 { 319 mp_small nv; 320 mp_result res = mp_int_to_int(n, &nv); 321 klisp_assert(res == MP_OK); 322 323 if (nv >= 0) 324 return mp_int_mul_pow2(K, x, nv, r); 325 else 326 return mp_int_div_pow2(K, x, -nv, r, NULL); 327 } 328 329 /* returns k, modifies all parameters (except f & p) */ 330 int32_t simple_fixup(klisp_State *K, Bigint *f, Bigint *p, Bigint *r, 331 Bigint *s, Bigint *mm, Bigint *mp) 332 { 333 mp_result res; 334 Bigint tmpz, tmpz2; 335 Bigint *tmp = &tmpz; 336 Bigint *tmp2 = &tmpz2; 337 Bigint onez; 338 Bigint *one = &onez; 339 res = mp_int_init(tmp); 340 res = mp_int_init(tmp2); 341 res = mp_int_init_value(K, one, 1); 342 res = mp_int_sub(K, p, one, tmp); 343 res = shift_2(K, one, tmp, tmp); 344 345 if (mp_int_compare(f, tmp) == 0) { 346 res = shift_2(K, mp, one, mp); 347 res = shift_2(K, r, one, r); 348 res = shift_2(K, s, one, s); 349 } 350 351 int k = 0; 352 353 /* tmp = ceiling (s/10), for while guard */ 354 res = mp_int_add_value(K, s, 9, tmp); 355 res = mp_int_div_value(K, tmp, 10, tmp, NULL); 356 357 while(mp_int_compare(r, tmp) < 0) { 358 --k; 359 res = mp_int_mul_value(K, r, 10, r); 360 res = mp_int_mul_value(K, mm, 10, mm); 361 res = mp_int_mul_value(K, mp, 10, mp); 362 /* tmp = ceiling (s/10), for while guard */ 363 res = mp_int_add_value(K, s, 9, tmp); 364 res = mp_int_div_value(K, tmp, 10, tmp, NULL); 365 } 366 367 /* tmp = 2r + mp; tmp2 = 2s */ 368 res = mp_int_mul_value(K, r, 2, tmp); 369 res = mp_int_add(K, tmp, mp, tmp); 370 res = mp_int_mul_value(K, s, 2, tmp2); 371 while(mp_int_compare(tmp, tmp2) >= 0) { 372 373 res = mp_int_mul_value(K, s, 10, s); 374 ++k; 375 376 /* tmp = 2r + mp; tmp2 = 2s */ 377 res = mp_int_mul_value(K, r, 2, tmp); 378 res = mp_int_add(K, tmp, mp, tmp); 379 res = mp_int_mul_value(K, s, 2, tmp2); 380 } 381 382 mp_int_clear(K, tmp); 383 mp_int_clear(K, tmp2); 384 mp_int_clear(K, one); 385 386 UNUSED(res); 387 return k; 388 } 389 390 /* TEMP: for now upoint is passed indicating where the least significant 391 integer digit should be (10^0 position) */ 392 #define digit_pos(k_, upoint_) ((k_) + (upoint_)) 393 394 bool dtoa(klisp_State *K, double d, char *buf, int32_t upoint, int32_t *out_h, 395 int32_t *out_k) 396 { 397 klisp_assert(sizeof(mp_small) == 4); 398 mp_result res; 399 Bigint e, p, f; 400 401 klisp_assert(d > 0.0); 402 403 /* convert d to three bigints m: significand, e: exponent & p: precision */ 404 /* d = m^(e-p) & m < 2^p */ 405 int ie, ip; 406 double mantissa = frexp(d, &ie); 407 ip = 0; 408 409 klisp_assert(mantissa * pow(2.0, ie) == d); 410 /* now 0.5 <= mantissa < 1 & mantissa * 2^expt = d */ 411 /* this could be a binary search, it could also be done reading the exponent 412 field of ieee754 directly... */ 413 while(mantissa != floor(mantissa)) { 414 mantissa *= 2.0; 415 ++ip; 416 } 417 418 /* mantissa is int & < 2^ip (was < 1=2^0 and by induction...) */ 419 klisp_assert(mantissa * pow(2.0, ie - ip) == d); 420 /* mantissa is at most 53 bits long as an int, load it in two parts 421 to f */ 422 int64_t im = (int64_t) mantissa; 423 /* f */ 424 /* cant load 32 bits at a time, second param is signed!, 425 but we know it's positive so load 32 then 31 */ 426 res = mp_int_init_value(K, &f, (mp_small) (im >> 31)); 427 res = mp_int_mul_pow2(K, &f, 31, &f); 428 res = mp_int_add_value(K, &f, (mp_small) im & 0x7fffffff, &f); 429 430 /* adjust f & p so that p is 53 TODO do in one step */ 431 /* XXX: is this is ok for denorms?? */ 432 while(ip < 53) { 433 ++ip; 434 res = mp_int_mul_value(K, &f, 2, &f); 435 } 436 437 /* e */ 438 res = mp_int_init_value(K, &e, (mp_small) ie); 439 440 /* p */ 441 res = mp_int_init_value(K, &p, (mp_small) ip); 442 443 /* start of FPP^2 algorithm */ 444 Bigint r, s; 445 Bigint mp, mm; 446 Bigint e_p; 447 Bigint zero, one; 448 449 res = mp_int_init_value(K, &zero, 0); 450 res = mp_int_init_value(K, &one, 1); 451 452 res = mp_int_init(&r); 453 res = mp_int_init(&s); 454 res = mp_int_init(&mm); 455 res = mp_int_init(&mp); 456 res = mp_int_init(&e_p); 457 458 res = mp_int_sub(K, &e, &p, &e_p); 459 460 // shift_2(f, max(e-p, 0), r); 461 // shift_2(1, max(-(e-p), 0), r); 462 if (mp_int_compare_zero(&e_p) >= 0) { 463 res = shift_2(K, &f, &e_p, &r); 464 res = shift_2(K, &one, &zero, &s); /* nop */ 465 res = shift_2(K, &one, &e_p, &mm); 466 } else { 467 res = shift_2(K, &f, &zero, &r); /* nop */ 468 res = mp_int_neg(K, &e_p, &e_p); 469 res = shift_2(K, &one, &e_p, &s); 470 res = shift_2(K, &one, &zero, &mm); 471 } 472 mp_int_copy(K, &mm, &mp); 473 474 int32_t k = simple_fixup(K, &f, &p, &r, &s, &mm, &mp); 475 int32_t h = k-1; 476 477 Bigint u, tmp, tmp2; 478 res = mp_int_init(&u); 479 res = mp_int_init(&tmp); 480 res = mp_int_init(&tmp2); 481 bool low, high; 482 483 do { 484 --k; 485 res = mp_int_mul_value(K, &r, 10, &tmp); 486 res = mp_int_div(K, &tmp, &s, &u, &r); 487 res = mp_int_mul_value(K, &mm, 10, &mm); 488 res = mp_int_mul_value(K, &mp, 10, &mp); 489 490 /* low/high flags */ 491 /* XXX try to make 1e23 round correctly, 492 it causes tmp == tmp2 but should probably 493 check oddness of digit and (may result in a digit 494 with value 10?, needing to backtrack) 495 In general make it so that if rounding done at reading 496 (should be round to even) is accounted for and the minimal 497 length number is generated */ 498 499 res = mp_int_mul_value(K, &r, 2, &tmp); 500 501 low = mp_int_compare(&tmp, &mm) < 0; 502 503 res = mp_int_mul_value(K, &s, 2, &tmp2); 504 res = mp_int_sub(K, &tmp2, &mp, &tmp2); 505 506 high = mp_int_compare(&tmp, &tmp2) > 0; 507 508 if (!low && !high) { 509 mp_small digit; 510 res = mp_int_to_int(&u, &digit); 511 klisp_assert(res == MP_OK); 512 klisp_assert(digit >= 0 && digit <= 9); 513 buf[digit_pos(k, upoint)] = '0' + digit; 514 } 515 } while(!low && !high); 516 517 mp_small digit; 518 res = mp_int_to_int(&u, &digit); 519 klisp_assert(res == MP_OK); 520 klisp_assert(digit >= 0 && digit <= 9); 521 522 if (low && high) { 523 res = mp_int_mul_value(K, &r, 2, &tmp); 524 int cmp = mp_int_compare(&tmp, &s); 525 if ((cmp == 0 && (digit & 1) != 0) || cmp > 0) 526 ++digit; 527 } else if (low) { 528 /* nothing */ 529 } else if (high) { 530 ++digit; 531 } else { 532 klisp_assert(0); 533 } 534 /* double check in case there was an increment */ 535 klisp_assert(digit >= 0 && digit <= 9); 536 buf[digit_pos(k, upoint)] = '0' + digit; 537 538 *out_h = h; 539 *out_k = k; 540 /* add '\0' to both sides */ 541 buf[digit_pos(k-1, upoint)] = '\0'; 542 buf[digit_pos(h+1, upoint)] = '\0'; 543 544 mp_int_clear(K, &f); 545 mp_int_clear(K, &e); 546 mp_int_clear(K, &p); 547 mp_int_clear(K, &r); 548 mp_int_clear(K, &s); 549 mp_int_clear(K, &mp); 550 mp_int_clear(K, &mm); 551 mp_int_clear(K, &e_p); 552 mp_int_clear(K, &zero); 553 mp_int_clear(K, &one); 554 mp_int_clear(K, &u); 555 mp_int_clear(K, &tmp); 556 mp_int_clear(K, &tmp2); 557 558 /* NOTE: digits are reversed! */ 559 return true; 560 } 561 562 563 /* TEMP: this is a stub for now, always return sufficiently large 564 number */ 565 int32_t kdouble_print_size(TValue tv_double) 566 { 567 klisp_assert(ttisdouble(tv_double)); 568 UNUSED(tv_double); 569 return 1024; 570 } 571 572 void kdouble_print_string(klisp_State *K, TValue tv_double, 573 char *buf, int32_t limit) 574 { 575 klisp_assert(ttisdouble(tv_double)); 576 /* TODO: add exponent to values too large or too small */ 577 /* TEMP */ 578 int32_t h = 0; 579 int32_t k = 0; 580 int32_t upoint = limit/2; 581 double od = dvalue(tv_double); 582 klisp_assert(!isnan(od) && !isinf(od)); 583 klisp_assert(limit > 10); 584 585 /* dtoa only works for d > 0 */ 586 if (od == 0.0) { 587 buf[0] = '0'; 588 buf[1] = '.'; 589 buf[2] = '0'; 590 buf[3] = '\0'; 591 return; 592 } 593 594 double d; 595 if (od < 0.0) 596 d = -od; 597 else d = od; 598 599 /* XXX this doesn't check limit, it should be large enough */ 600 UNUSED(dtoa(K, d, buf, upoint, &h, &k)); 601 602 klisp_assert(upoint + k >= 0 && upoint + h + 1 < limit); 603 604 /* rearrange the digits */ 605 /* TODO do this more efficiently */ 606 607 608 int32_t size = h - k + 1; 609 int32_t start = upoint+k; 610 /* first reverse the digits */ 611 for (int32_t i = upoint+k, j = upoint+h; i < j; i++, j--) { 612 char ch = buf[i]; 613 buf[i] = buf[j]; 614 buf[j] = ch; 615 } 616 617 /* TODO use exponents */ 618 619 /* if necessary make room for leading zeros and sign, 620 move all to the left */ 621 622 int extra_size = (od < 0? 1 : 0) + (h < 0? 2 + (-h-1) : 0); 623 624 klisp_assert(extra_size + size + 2 < limit); 625 memmove(buf+extra_size, buf+start, size); 626 627 int32_t i = 0; 628 if (od < 0) 629 buf[i++] = '-'; 630 631 if (h < 0) { 632 /* fraction with leading 0. and with possibly more leading zeros */ 633 buf[i++] = '0'; 634 buf[i++] = '.'; 635 for (int32_t j = -1; j > h; j--) { 636 buf[i++] = '0'; 637 } 638 int frac_size = size; 639 i += frac_size; 640 buf[i++] = '\0'; 641 } else if (k >= 0) { 642 /* integer with possibly trailing zeros */ 643 klisp_assert(size+extra_size+k+4 < limit); 644 int int_size = size; 645 i += int_size; 646 for (int32_t j = 0; j < k; j++) { 647 buf[i++] = '0'; 648 } 649 buf[i++] = '.'; 650 buf[i++] = '0'; 651 buf[i++] = '\0'; 652 } else { /* both integer and fractional part, make room for the point */ 653 /* k < 0, h >= 0 */ 654 int32_t int_size = h+1; 655 int32_t frac_size = -k; 656 memmove(buf+i+int_size+1, buf+i+int_size, frac_size); 657 i += int_size; 658 buf[i++] = '.'; 659 i += frac_size; 660 buf[i++] = '\0'; 661 } 662 return; 663 } 664 665 double kdouble_div_mod(double n, double d, double *res_mod) 666 { 667 double div = floor(n / d); 668 double mod = fmod(n, d); 669 670 /* div, mod or div-and-mod */ 671 /* 0 <= mod0 < |d| */ 672 if (mod < 0.0) { 673 if (d < 0.0) { 674 mod -= d; 675 div += 1.0; 676 } else { 677 mod += d; 678 div -= 1.0; 679 } 680 } 681 *res_mod = mod; 682 return div; 683 } 684 685 double kdouble_div0_mod0(double n, double d, double *res_mod) 686 { 687 double div = floor(n / d); 688 double mod = fmod(n, d); 689 690 /* div0, mod0 or div-and-mod0 */ 691 /* 692 ** Adjust q and r so that: 693 ** -|d/2| <= mod0 < |d/2| which is the same as 694 ** dmin <= mod0 < dmax, where 695 ** dmin = -|d/2| and dmax = |d/2| 696 */ 697 double dmin = -((d<0.0? -d : d) / 2.0); 698 double dmax = (d<0.0? -d : d) / 2.0; 699 700 if (mod < dmin) { 701 if (d < 0) { 702 mod -= d; 703 div += 1.0; 704 } else { 705 mod += d; 706 div -= 1.0; 707 } 708 } else if (mod >= dmax) { 709 if (d < 0) { 710 mod += d; 711 div += 1.0; 712 } else { 713 mod -= d; 714 div -= 1.0; 715 } 716 } 717 *res_mod = mod; 718 return div; 719 } 720 721 TValue kdouble_to_integer(klisp_State *K, TValue tv_double, kround_mode mode) 722 { 723 double d = dvalue(tv_double); 724 switch(mode) { 725 case K_TRUNCATE: 726 d = trunc(d); 727 break; 728 case K_CEILING: 729 d = ceil(d); 730 break; 731 case K_FLOOR: 732 d = floor(d); 733 break; 734 case K_ROUND_EVEN: { 735 int res = fesetround(FE_TONEAREST); /* REFACTOR: should be done once only... */ 736 klisp_assert(res == 0); 737 d = nearbyint(d); 738 break; 739 } 740 } 741 /* ASK John: we currently return inexact if given inexact is this ok? 742 or should it return an exact integer? */ 743 return ktag_double(d); 744 #if 0 745 tv_double = ktag_double(d); /* won't alloc mem so no need to root */ 746 return kinexact_to_exact(K, tv_double); 747 #endif 748 }