source: clfs-embedded/patches/mpfr-2.3.2-branch_update-1.patch@ 128e9a4

Last change on this file since 128e9a4 was 06beb7a, checked in by Jim Gifford <clfs@…>, 16 years ago

Added: Patches
  • Property mode set to 100644
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  
    3131#define MPFR_MAX_BASE 62
    3232
    3333struct 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) */
    3945};
    4046
    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*/
    4252static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = {
    4353  {1UL, 1UL},
    4454  {53UL, 84UL},
    … …  
    441451  MPFR_ZIV_DECL (loop);
    442452  MPFR_TMP_DECL (marker);
    443453
    444   /* determine the minimal precision for the computation */
     454  /* initialize the working precision */
    445455  prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x));
    446456
    447457  /* compute y as long as rounding is not possible */
    … …  
    449459  MPFR_ZIV_INIT (loop, prec);
    450460  for (;;)
    451461    {
    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). */
    453465      ysize = (prec - 1) / BITS_PER_MP_LIMB + 1;
     466      /* prec bits corresponds to ysize limbs */
    454467      ysize_bits = ysize * BITS_PER_MP_LIMB;
     468      /* and to ysize_bits >= prec > MPFR_PREC (x) bits */
    455469      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 */
    457471
    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 =>
    462478          ysize*BITS_PER_MP_LIMB can not overflow.
    463          Compute pstr_size = ysize_bits * Num / Den
    464           Where Num/Den ~ 1/log2(base)
     479         We compute pstr_size = 1 + ceil(ysize_bits * Num / Den)
     480          where Num/Den >= 1/log2(base)
    465481         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) */
    473500
     501      /* if the number of wanted characters is more than what we have in
     502         pstr->mant, round it down */
    474503      if (pstr_size >= pstr->prec)
    475504        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);
    477506
    478       /* convert str into binary */
     507      /* convert str into binary: note that pstr->mant is big endian,
     508         thus no offset is needed */
    479509      real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base);
    480       MPFR_ASSERTD ( real_ysize <= ysize+1);
     510      MPFR_ASSERTD (real_ysize <= ysize+1);
    481511
    482512      /* 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 */
    484514      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 */
    487519        {
     520          /* we have enough limbs to store {y, real_ysize} */
    488521          /* shift {y, num_limb} for count bits to the left */
    489522          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);
    492524          if (real_ysize != ysize)
    493525            {
    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);
    495528              MPN_ZERO (y, ysize - real_ysize);
    496529            }
    497530          /* for each bit shift decrease exponent of y */
    498531          /* (This should not overflow) */
    499532          exp = - ((ysize - real_ysize) * BITS_PER_MP_LIMB + count);
    500533        }
    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) */
    502537        {
    503538          /* 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;
    506544          /* for each bit shift increase exponent of y */
    507545          exp = BITS_PER_MP_LIMB - count;
    508546        }
    … …  
    542580          result = y;
    543581          err = 0;
    544582        }
    545       /* case pstr->exp_base > pstr_size */
     583      /* case non-power-of-two-base, and pstr->exp_base > pstr_size */
    546584      else if (pstr->exp_base > (mp_exp_t) pstr_size)
    547585        {
    548586          mp_limb_t *z;
    … …  
    559597            goto overflow;
    560598          exact = exact && (err == -1);
    561599
     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
    562607          /* multiply(y = 0.mant_s[0]...mant_s[pr-1])_base by base^(exp_s-g) */
    563608          mpn_mul_n (result, y, z, ysize);
    564609
    … …  
    588633              exp --;
    589634            }
    590635
     636          /* if the low ysize limbs of {result, 2*ysize} are all zero,
     637             then the result is still "exact" (if it was before) */
    591638          exact = exact && (mpn_scan1 (result, 0)
    592639                            >= (unsigned long) ysize_bits);
    593640          result += ysize;
    … …  
    598645          mp_limb_t *z;
    599646          mp_exp_t exp_z;
    600647
    601           result = (mp_limb_t*) MPFR_TMP_ALLOC ( (3*ysize+1) * BYTES_PER_MP_LIMB);
     648          result = (mp_limb_t*) MPFR_TMP_ALLOC ((3*ysize+1) * BYTES_PER_MP_LIMB);
    602649
    603650          /* set y to y * K^ysize */
    604651          y = y - ysize;  /* we have allocated ysize limbs at y - ysize */
    … …  
    622669          /* compute y / z */
    623670          /* result will be put into result + n, and remainder into result */
    624671          mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y,
    625                        2*ysize, z, ysize);
     672                       2 * ysize, z, ysize);
    626673
    627674          /* exp -= exp_z + ysize_bits with overflow checking
    628675             and check that we can add/substract 2 to exp without overflow */
    … …  
    635682                              MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
    636683                              goto overflow, goto underflow);
    637684          err += 2;
     685          /* if the remainder of the division is zero, then the result is
     686             still "exact" if it was before */
    638687          exact = exact && (mpn_popcount (result, ysize) == 0);
    639688
    640689          /* normalize result */
    … …  
    648697            }
    649698          result += ysize;
    650699        }
    651       /* case exp_base = pstr_size */
     700      /* case exp_base = pstr_size: no multiplication or division needed */
    652701      else
    653702        {
    654703          /* base^(exp_s-pr) = 1             nothing to compute */
    … …  
    656705          err = 0;
    657706        }
    658707
     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
    659719      /* test if rounding is possible, and if so exit the loop */
    660720      if (exact || mpfr_can_round_raw (result, ysize,
    661721                                       (pstr->negative) ? -1 : 1,
    … …  
    679739      exp ++;
    680740    }
    681741
     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
    682749  /* Set sign of x before exp since check_range needs a valid sign */
    683750  (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
    684751

  • 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  
    8181  mpfr_clear (a);
    8282}
    8383
     84/* Bug found by Christoph Lauter. */
     85static void
     86bug20081028 (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
    84102int
    85103main (int argc, char *argv[])
    86104{
    … …  
    845863  mpfr_clear (y);
    846864
    847865  check_underflow ();
     866  bug20081028 ();
    848867
    849868  tests_end_mpfr ();
    850869  return 0;

  • tests/tstrtofr.c 

    diff -Naur mpfr-2.3.2.orig/tests/tstrtofr.c mpfr-2.3.2/tests/tstrtofr.c
    old new  
    2727
    2828#include "mpfr-test.h"
    2929
    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
     30static void
    5231check_special (void)
    5332{
    5433  mpfr_t x, y;
    … …  
    551530"1.001000010110011011000101100000101111101001100011101101001111110111000010010110010001100e-16920"}
    552531};
    553532
    554 
    555 void
     533static void
    556534check_reftable ()
    557535{
    558536  int i, base;
    … …  
    597575  mpfr_clear (x);
    598576}
    599577
    600 void
     578static void
    601579check_parse (void)
    602580{
    603581  mpfr_t x;
    … …  
    951929
    952930  mpfr_clear (x);
    953931}
     932
     933/* Bug found by Christoph Lauter (in mpfr_set_str). */
     934static 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
     968static void
     969bug20081028 (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
     1018int
     1019main (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.