source:
clfs-embedded/patches/mpfr-2.3.2-branch_update-1.patch@
88b7bb1
Last change on this file since 88b7bb1 was 06beb7a, checked in by , 16 years ago | |
---|---|
|
|
File size: 16.4 KB |
-
strtofr.c
Submitted By: Jim Gifford (jim at cross-lfs dot org) Date: 2009-01-05 Initial Package Version: 2.3.2 Origin: MPFR Website Upstream Status: Fixed Description: See http://www.mpfr.org Website Under Bugs diff -Naur mpfr-2.3.2.orig/strtofr.c mpfr-2.3.2/strtofr.c
old new 31 31 #define MPFR_MAX_BASE 62 32 32 33 33 struct parsed_string { 34 int negative; 35 int base; 36 unsigned char *mantissa, *mant; 37 size_t prec, alloc; 38 mp_exp_t exp_base, exp_bin; 34 int negative; /* non-zero iff the number is negative */ 35 int base; /* base of the string */ 36 unsigned char *mantissa; /* raw significand (without any point) */ 37 unsigned char *mant; /* stripped significand (without starting and 38 ending zeroes) */ 39 size_t prec; /* length of mant (zero for +/-0) */ 40 size_t alloc; /* allocation size of mantissa */ 41 mp_exp_t exp_base; /* number of digits before the point */ 42 mp_exp_t exp_bin; /* exponent in case base=2 or 16, and the pxxx 43 format is used (i.e., exponent is given in 44 base 10) */ 39 45 }; 40 46 41 /* This table has been generated by the following program */ 47 /* This table has been generated by the following program. 48 For 2 <= b <= MPFR_MAX_BASE, 49 RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1] 50 is an upper approximation of log(2)/log(b). 51 */ 42 52 static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = { 43 53 {1UL, 1UL}, 44 54 {53UL, 84UL}, … … 441 451 MPFR_ZIV_DECL (loop); 442 452 MPFR_TMP_DECL (marker); 443 453 444 /* determine the minimal precision for the computation */454 /* initialize the working precision */ 445 455 prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); 446 456 447 457 /* compute y as long as rounding is not possible */ … … 449 459 MPFR_ZIV_INIT (loop, prec); 450 460 for (;;) 451 461 { 452 /* initialize y to the value of 0.mant_s[0]...mant_s[pr-1] */ 462 /* Set y to the value of the ~prec most significant bits of pstr->mant 463 (as long as we guarantee correct rounding, we don't need to get 464 exactly prec bits). */ 453 465 ysize = (prec - 1) / BITS_PER_MP_LIMB + 1; 466 /* prec bits corresponds to ysize limbs */ 454 467 ysize_bits = ysize * BITS_PER_MP_LIMB; 468 /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ 455 469 y = (mp_limb_t*) MPFR_TMP_ALLOC ((2 * ysize + 1) * sizeof (mp_limb_t)); 456 y += ysize; 470 y += ysize; /* y has (ysize+1) allocated limbs */ 457 471 458 /* required precision for str to have ~ysize limbs 459 We must have (2^(BITS_PER_MP_LIMB))^ysize ~= base^pstr_size 460 aka pstr_size = ceil (ysize*BITS_PER_MP_LIMB/log2(base)) 461 ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => 472 /* pstr_size is the number of characters we read in pstr->mant 473 to have at least ysize full limbs. 474 We must have base^(pstr_size-1) >= (2^(BITS_PER_MP_LIMB))^ysize 475 (in the worst case, the first digit is one and all others are zero). 476 i.e., pstr_size >= 1 + ysize*BITS_PER_MP_LIMB/log2(base) 477 Since ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => 462 478 ysize*BITS_PER_MP_LIMB can not overflow. 463 Compute pstr_size = ysize_bits * Num / Den464 Where Num/Den ~1/log2(base)479 We compute pstr_size = 1 + ceil(ysize_bits * Num / Den) 480 where Num/Den >= 1/log2(base) 465 481 It is not exactly ceil(1/log2(base)) but could be one more (base 2) 466 Quite ugly since it tries to avoid overflow */ 467 pstr_size = ( ysize_bits / RedInvLog2Table[pstr->base-2][1] 468 * RedInvLog2Table[pstr->base-2][0] ) 469 + (( ysize_bits % RedInvLog2Table[pstr->base-2][1]) 470 * RedInvLog2Table[pstr->base-2][0] 471 / RedInvLog2Table[pstr->base-2][1] ) 472 + 1; 482 Quite ugly since it tries to avoid overflow: 483 let Num = RedInvLog2Table[pstr->base-2][0] 484 and Den = RedInvLog2Table[pstr->base-2][1], 485 and ysize_bits = a*Den+b, 486 then ysize_bits * Num/Den = a*Num + (b * Num)/Den, 487 thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den 488 */ 489 { 490 unsigned long Num = RedInvLog2Table[pstr->base-2][0]; 491 unsigned long Den = RedInvLog2Table[pstr->base-2][1]; 492 pstr_size = ((ysize_bits / Den) * Num) 493 + (((ysize_bits % Den) * Num + Den - 1) / Den) 494 + 1; 495 } 496 497 /* since pstr_size corresponds to at least ysize_bits full bits, 498 and ysize_bits > prec, the weight of the neglected part of 499 pstr->mant (if any) is < ulp(y) < ulp(x) */ 473 500 501 /* if the number of wanted characters is more than what we have in 502 pstr->mant, round it down */ 474 503 if (pstr_size >= pstr->prec) 475 504 pstr_size = pstr->prec; 476 MPFR_ASSERTD ( (mp_exp_t)pstr_size == (mp_exp_t) pstr_size);505 MPFR_ASSERTD (pstr_size == (mp_exp_t) pstr_size); 477 506 478 /* convert str into binary */ 507 /* convert str into binary: note that pstr->mant is big endian, 508 thus no offset is needed */ 479 509 real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base); 480 MPFR_ASSERTD ( 510 MPFR_ASSERTD (real_ysize <= ysize+1); 481 511 482 512 /* normalize y: warning we can get even get ysize+1 limbs! */ 483 MPFR_ASSERTD (y[real_ysize - 1] != 0); 513 MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ 484 514 count_leading_zeros (count, y[real_ysize - 1]); 485 exact = (real_ysize <= ysize); 486 if (exact != 0) /* shift y to the left in that case y should be exact */ 515 /* exact means that the number of limbs of the output of mpn_set_str 516 is less or equal to ysize */ 517 exact = real_ysize <= ysize; 518 if (exact) /* shift y to the left in that case y should be exact */ 487 519 { 520 /* we have enough limbs to store {y, real_ysize} */ 488 521 /* shift {y, num_limb} for count bits to the left */ 489 522 if (count != 0) 490 mpn_lshift (y, y, real_ysize, count); 491 /* shift {y, num_limb} for (ysize-num_limb) limbs to the left */ 523 mpn_lshift (y + ysize - real_ysize, y, real_ysize, count); 492 524 if (real_ysize != ysize) 493 525 { 494 MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); 526 if (count == 0) 527 MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); 495 528 MPN_ZERO (y, ysize - real_ysize); 496 529 } 497 530 /* for each bit shift decrease exponent of y */ 498 531 /* (This should not overflow) */ 499 532 exp = - ((ysize - real_ysize) * BITS_PER_MP_LIMB + count); 500 533 } 501 else /* shift y for the right */ 534 else /* shift y to the right, by doing this we might lose some 535 bits from the result of mpn_set_str (in addition to the 536 characters neglected from pstr->mant) */ 502 537 { 503 538 /* shift {y, num_limb} for (BITS_PER_MP_LIMB - count) bits 504 to the right */ 505 mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count); 539 to the right. FIXME: can we prove that count cannot be zero here, 540 since mpn_rshift does not accept a shift of BITS_PER_MP_LIMB? */ 541 MPFR_ASSERTD (count != 0); 542 exact = mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count) == 543 MPFR_LIMB_ZERO; 506 544 /* for each bit shift increase exponent of y */ 507 545 exp = BITS_PER_MP_LIMB - count; 508 546 } … … 542 580 result = y; 543 581 err = 0; 544 582 } 545 /* case pstr->exp_base > pstr_size */583 /* case non-power-of-two-base, and pstr->exp_base > pstr_size */ 546 584 else if (pstr->exp_base > (mp_exp_t) pstr_size) 547 585 { 548 586 mp_limb_t *z; … … 559 597 goto overflow; 560 598 exact = exact && (err == -1); 561 599 600 /* If exact is non zero, then z equals exactly the value of the 601 pstr_size most significant digits from pstr->mant, i.e., the 602 only difference can come from the neglected pstr->prec-pstr_size 603 least significant digits of pstr->mant. 604 If exact is zero, then z is rounded towards zero with respect 605 to that value. */ 606 562 607 /* multiply(y = 0.mant_s[0]...mant_s[pr-1])_base by base^(exp_s-g) */ 563 608 mpn_mul_n (result, y, z, ysize); 564 609 … … 588 633 exp --; 589 634 } 590 635 636 /* if the low ysize limbs of {result, 2*ysize} are all zero, 637 then the result is still "exact" (if it was before) */ 591 638 exact = exact && (mpn_scan1 (result, 0) 592 639 >= (unsigned long) ysize_bits); 593 640 result += ysize; … … 598 645 mp_limb_t *z; 599 646 mp_exp_t exp_z; 600 647 601 result = (mp_limb_t*) MPFR_TMP_ALLOC ( 648 result = (mp_limb_t*) MPFR_TMP_ALLOC ((3*ysize+1) * BYTES_PER_MP_LIMB); 602 649 603 650 /* set y to y * K^ysize */ 604 651 y = y - ysize; /* we have allocated ysize limbs at y - ysize */ … … 622 669 /* compute y / z */ 623 670 /* result will be put into result + n, and remainder into result */ 624 671 mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y, 625 2 *ysize, z, ysize);672 2 * ysize, z, ysize); 626 673 627 674 /* exp -= exp_z + ysize_bits with overflow checking 628 675 and check that we can add/substract 2 to exp without overflow */ … … 635 682 MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, 636 683 goto overflow, goto underflow); 637 684 err += 2; 685 /* if the remainder of the division is zero, then the result is 686 still "exact" if it was before */ 638 687 exact = exact && (mpn_popcount (result, ysize) == 0); 639 688 640 689 /* normalize result */ … … 648 697 } 649 698 result += ysize; 650 699 } 651 /* case exp_base = pstr_size */700 /* case exp_base = pstr_size: no multiplication or division needed */ 652 701 else 653 702 { 654 703 /* base^(exp_s-pr) = 1 nothing to compute */ … … 656 705 err = 0; 657 706 } 658 707 708 /* If result is exact, we still have to consider the neglected part 709 of the input string. For a directed rounding, in that case we could 710 still correctly round, since the neglected part is less than 711 one ulp, but that would make the code more complex, and give a 712 speedup for rare cases only. */ 713 exact = exact && (pstr_size == pstr->prec); 714 715 /* at this point, result is an approximation rounded towards zero 716 of the pstr_size most significant digits of pstr->mant, with 717 equality in case exact is non-zero. */ 718 659 719 /* test if rounding is possible, and if so exit the loop */ 660 720 if (exact || mpfr_can_round_raw (result, ysize, 661 721 (pstr->negative) ? -1 : 1, … … 679 739 exp ++; 680 740 } 681 741 742 if (res == 0) /* fix ternary value */ 743 { 744 exact = exact && (pstr_size == pstr->prec); 745 if (!exact) 746 res = (pstr->negative) ? 1 : -1; 747 } 748 682 749 /* Set sign of x before exp since check_range needs a valid sign */ 683 750 (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); 684 751 -
tests/tset_str.c
diff -Naur mpfr-2.3.2.orig/tests/tset_str.c mpfr-2.3.2/tests/tset_str.c
old new 81 81 mpfr_clear (a); 82 82 } 83 83 84 /* Bug found by Christoph Lauter. */ 85 static void 86 bug20081028 (void) 87 { 88 mpfr_t x; 89 const char *s = "0.10000000000000000000000000000001E1"; 90 91 mpfr_init2 (x, 32); 92 mpfr_set_str (x, "1.00000000000000000006", 10, GMP_RNDU); 93 if (! mpfr_greater_p (x, __gmpfr_one)) 94 { 95 printf ("Error in bug20081028:\nExpected %s\nGot ", s); 96 mpfr_dump (x); 97 exit (1); 98 } 99 mpfr_clear (x); 100 } 101 84 102 int 85 103 main (int argc, char *argv[]) 86 104 { … … 845 863 mpfr_clear (y); 846 864 847 865 check_underflow (); 866 bug20081028 (); 848 867 849 868 tests_end_mpfr (); 850 869 return 0; -
tests/tstrtofr.c
diff -Naur mpfr-2.3.2.orig/tests/tstrtofr.c mpfr-2.3.2/tests/tstrtofr.c
old new 27 27 28 28 #include "mpfr-test.h" 29 29 30 static void check_reftable (void); 31 static void check_special (void); 32 static void check_retval (void); 33 static void check_overflow (void); 34 static void check_parse (void); 35 36 int 37 main (int argc, char *argv[]) 38 { 39 tests_start_mpfr (); 40 41 check_special(); 42 check_reftable (); 43 check_parse (); 44 check_overflow (); 45 check_retval (); 46 47 tests_end_mpfr (); 48 return 0; 49 } 50 51 void 30 static void 52 31 check_special (void) 53 32 { 54 33 mpfr_t x, y; … … 551 530 "1.001000010110011011000101100000101111101001100011101101001111110111000010010110010001100e-16920"} 552 531 }; 553 532 554 555 void 533 static void 556 534 check_reftable () 557 535 { 558 536 int i, base; … … 597 575 mpfr_clear (x); 598 576 } 599 577 600 void578 static void 601 579 check_parse (void) 602 580 { 603 581 mpfr_t x; … … 951 929 952 930 mpfr_clear (x); 953 931 } 932 933 /* Bug found by Christoph Lauter (in mpfr_set_str). */ 934 static struct bug20081025_test { 935 mpfr_rnd_t rnd; 936 int inexact; 937 const char *str; 938 const char *binstr; 939 } Bug20081028Table[] = { 940 {GMP_RNDN, -1, "1.00000000000000000006", "1"}, 941 {GMP_RNDZ, -1, "1.00000000000000000006", "1"}, 942 {GMP_RNDU, +1, "1.00000000000000000006", 943 "10000000000000000000000000000001e-31"}, 944 {GMP_RNDD, -1, "1.00000000000000000006", "1"}, 945 946 947 {GMP_RNDN, +1, "-1.00000000000000000006", "-1"}, 948 {GMP_RNDZ, +1, "-1.00000000000000000006", "-1"}, 949 {GMP_RNDU, +1, "-1.00000000000000000006", "-1"}, 950 {GMP_RNDD, -1, "-1.00000000000000000006", 951 "-10000000000000000000000000000001e-31"}, 952 953 {GMP_RNDN, +1, "0.999999999999999999999999999999999999999999999", "1"}, 954 {GMP_RNDZ, -1, "0.999999999999999999999999999999999999999999999", 955 "11111111111111111111111111111111e-32"}, 956 {GMP_RNDU, +1, "0.999999999999999999999999999999999999999999999", "1"}, 957 {GMP_RNDD, -1, "0.999999999999999999999999999999999999999999999", 958 "11111111111111111111111111111111e-32"}, 959 960 {GMP_RNDN, -1, "-0.999999999999999999999999999999999999999999999", "-1"}, 961 {GMP_RNDZ, +1, "-0.999999999999999999999999999999999999999999999", 962 "-11111111111111111111111111111111e-32"}, 963 {GMP_RNDU, +1, "-0.999999999999999999999999999999999999999999999", 964 "-11111111111111111111111111111111e-32"}, 965 {GMP_RNDD, -1, "-0.999999999999999999999999999999999999999999999", "-1"} 966 }; 967 968 static void 969 bug20081028 (void) 970 { 971 int i; 972 int inexact, res; 973 mpfr_rnd_t rnd; 974 mpfr_t x, y; 975 char *s; 976 977 mpfr_init2 (x, 32); 978 mpfr_init2 (y, 32); 979 for (i = 0 ; i < numberof (Bug20081028Table) ; i++) 980 { 981 rnd = Bug20081028Table[i].rnd; 982 inexact = Bug20081028Table[i].inexact; 983 mpfr_set_str_binary (x, Bug20081028Table[i].binstr); 984 res = mpfr_strtofr (y, Bug20081028Table[i].str, &s, 10, rnd); 985 if (s == NULL || *s != 0) 986 { 987 printf ("Error in Bug20081028: strtofr didn't parse entire input\n" 988 "for (i=%d) Str=\"%s\"", i, Bug20081028Table[i].str); 989 exit (1); 990 } 991 if (! SAME_SIGN (res, inexact)) 992 { 993 printf ("Error in Bug20081028: expected %s ternary value, " 994 "got %d\nfor (i=%d) Rnd=%s Str=\"%s\"\n Set binary gives: ", 995 inexact > 0 ? "positive" : "negative", 996 res, i, mpfr_print_rnd_mode(rnd), Bug20081028Table[i].str); 997 mpfr_dump (x); 998 printf (" strtofr gives: "); 999 mpfr_dump (y); 1000 exit (1); 1001 } 1002 if (mpfr_cmp (x, y)) 1003 { 1004 printf ("Error in Bug20081028: Results differ between strtofr and " 1005 "set_binary\nfor (i=%d) Rnd=%s Str=\"%s\"\n" 1006 " Set binary gives: ", 1007 i, mpfr_print_rnd_mode(rnd), Bug20081028Table[i].str); 1008 mpfr_dump (x); 1009 printf (" strtofr gives: "); 1010 mpfr_dump (y); 1011 exit (1); 1012 } 1013 } 1014 mpfr_clear (y); 1015 mpfr_clear (x); 1016 } 1017 1018 int 1019 main (int argc, char *argv[]) 1020 { 1021 tests_start_mpfr (); 1022 1023 check_special(); 1024 check_reftable (); 1025 check_parse (); 1026 check_overflow (); 1027 check_retval (); 1028 bug20081028 (); 1029 1030 tests_end_mpfr (); 1031 return 0; 1032 }
Note:
See TracBrowser
for help on using the repository browser.