Ignore:
Timestamp:
Jan 10, 2006, 8:03:34 AM (13 years ago)
Author:
mww
Message:

Bug: #6237
Submitted by: vincent-opdarw@…
Reviewed by: mww@

use latest patch-set, inc. revision

Location:
trunk/dports/devel/mpfr
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/dports/devel/mpfr/Portfile

    r14408 r15794  
    1 # $Id: Portfile,v 1.9 2005/10/04 12:31:34 mww Exp $
     1# $Id: Portfile,v 1.10 2006/01/10 08:03:33 mww Exp $
    22
    33PortSystem 1.0
     
    55name            mpfr
    66version         2.2.0
    7 revision        2
     7revision        3
    88categories      devel math
    99platforms       darwin
  • trunk/dports/devel/mpfr/files/patches

    r14408 r15794  
    1 diff -Naurd lngamma.c lngamma.c
     1diff -Naurd mpfr-2.2.0/lngamma.c mpfr-2.2.0-p1/lngamma.c
    22--- lngamma.c   2005-09-09 15:17:58.000000000 +0000
    33+++ lngamma.c   2005-09-29 11:27:04.000000000 +0000
     
    1313 
    1414   mpfr_init2 (s, MPFR_PREC_MIN);
    15 diff -Naurd tests/tlngamma.c tests/tlngamma.c
     15diff -Naurd mpfr-2.2.0/tests/tlngamma.c mpfr-2.2.0-p1/tests/tlngamma.c
    1616--- tests/tlngamma.c    2005-09-09 15:17:59.000000000 +0000
    1717+++ tests/tlngamma.c    2005-09-29 11:20:34.000000000 +0000
     
    4040   mpfr_set_prec (y, 53);
    4141 
    42 diff -Naurd mpfr.h mpfr.h
     42diff -Naurd mpfr-2.2.0-p1/mpfr.h mpfr-2.2.0-p2/mpfr.h
    4343--- mpfr.h      2005-09-06 15:02:12.000000000 +0000
    4444+++ mpfr.h      2005-09-29 11:36:36.000000000 +0000
     
    6464 #define mpfr_set_si(_f,_s,_r)              \
    6565  (__builtin_constant_p (_s) && (_s) >= 0 ? \
    66 diff -Naurd tests/tset_si.c tests/tset_si.c
     66diff -Naurd mpfr-2.2.0-p1/tests/tset_si.c mpfr-2.2.0-p2/tests/tset_si.c
    6767--- tests/tset_si.c     2005-08-18 17:03:17.000000000 +0000
    6868+++ tests/tset_si.c     2005-09-29 09:19:39.000000000 +0000
     
    111111   return 0;
    112112 }
    113 diff -Naurd configure configure
     113diff -Naurd mpfr-2.2.0-p2/configure mpfr-2.2.0-p3/configure
    114114--- configure   2005-09-19 13:31:58.000000000 +0000
    115115+++ configure   2005-10-02 10:49:55.000000000 +0000
     
    123123       case "$host_os" in
    124124       rhapsody* | darwin1.[012])
     125diff -Naurd mpfr-2.2.0-p3/tests/tpow.c mpfr-2.2.0-p4/tests/tpow.c
     126--- tests/tpow.c        2005-06-02 16:12:05.000000000 +0000
     127+++ tests/tpow.c        2005-10-06 09:54:52.000000000 +0000
     128@@ -509,6 +509,7 @@
     129   for (i = 0; i < 11; i++)
     130     for (j = 0; j < 11; j++)
     131       {
     132+        double d;
     133         int p;
     134         static int q[11][11] = {
     135           /*          NaN +inf -inf  +0   -0   +1   -1   +2   -2  +0.5 -0.5 */
     136@@ -527,7 +528,7 @@
     137         test_pow (r, t[i], t[j], GMP_RNDN);
     138         p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
     139           mpfr_cmp_ui (r, 0) == 0 ? 2 :
     140-          (int) (fabs (mpfr_get_d (r, GMP_RNDN)) * 128.0);
     141+          (d = mpfr_get_d (r, GMP_RNDN), (int) (ABS(d) * 128.0));
     142         if (p != 0 && MPFR_SIGN(r) < 0)
     143           p = -p;
     144         if (p != q[i][j])
     145diff -Naurd mpfr-2.2.0-p4/cache.c mpfr-2.2.0-p5/cache.c
     146--- cache.c     2005-08-18 16:35:03.000000000 +0000
     147+++ cache.c     2005-11-23 09:04:29.000000000 +0000
     148@@ -32,9 +32,9 @@
     149 void
     150 mpfr_clear_cache (mpfr_cache_t cache)
     151 {
     152-  if (MPFR_PREC(cache->x) != 0)
     153+  if (MPFR_PREC (cache->x) != 0)
     154     mpfr_clear (cache->x);
     155-  MPFR_PREC(cache->x) = 0;
     156+  MPFR_PREC (cache->x) = 0;
     157 }
     158 
     159 int
     160@@ -47,45 +47,58 @@
     161 
     162   MPFR_SAVE_EXPO_MARK (expo);
     163 
     164-  /* Check if the cache has been already filled */
     165-  if (MPFR_UNLIKELY(pold == 0))
     166-    mpfr_init2 (cache->x, prec);
     167+  if (MPFR_UNLIKELY (prec > pold))
     168+    {
     169+      /* No previous result in the cache or the precision of the
     170+         previous result is not sufficient. */
     171 
     172-  /* Check if we can round with the previous result */
     173-  else if (MPFR_LIKELY(prec <= pold))
     174-    goto round;
     175+      if (MPFR_UNLIKELY (pold == 0))  /* No previous result. */
     176+        mpfr_init2 (cache->x, prec);
     177 
     178-  /* Update the cache */
     179-  pold = prec /*MPFR_PREC_MIN + prec + __gmpfr_ceil_exp2 (prec)*/;
     180-  MPFR_ASSERTD (pold >= prec);
     181-  mpfr_prec_round (cache->x, pold, GMP_RNDN);
     182-  cache->inexact = (*cache->func) (cache->x, GMP_RNDN);
     183+      /* Update the cache. */
     184+      pold = prec;
     185+      mpfr_prec_round (cache->x, pold, GMP_RNDN);
     186+      cache->inexact = (*cache->func) (cache->x, GMP_RNDN);
     187+    }
     188 
     189- round:
     190-  /* First check if the cache has the exact value (Unlikely)
     191-     Else the exact value is between (assuming x=cache->x > 0)
     192-     x and x+ulp(x) if cache->inexact < 0
     193-     x-ulp(x) and x if cache->inexact > 0
     194-     and abs(x-exact) <= ulp(x)/2 */
     195-  MPFR_ASSERTD (MPFR_IS_POS(cache->x)); /* TODO...*/
     196-  /* We must use nextbelow instead of sub_one_ulp, since we know
     197-     that the exact value is < 1/2ulp(x) (We want sub_demi_ulp(x)).
     198-     Can't use mpfr_set since we need the even flag. */
     199+  /* First, check if the cache has the exact value (unlikely).
     200+     Else the exact value is between (assuming x=cache->x > 0):
     201+       x and x+ulp(x) if cache->inexact < 0,
     202+       x-ulp(x) and x if cache->inexact > 0,
     203+     and abs(x-exact) <= ulp(x)/2. */
     204+  MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */
     205   sign = MPFR_SIGN (cache->x);
     206   MPFR_SET_EXP (dest, MPFR_GET_EXP (cache->x));
     207   MPFR_SET_SIGN (dest, sign);
     208-  MPFR_RNDRAW_EVEN (inexact, dest,
     209-                    MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign,
     210-                    if (MPFR_UNLIKELY ( ++MPFR_EXP (dest) > __gmpfr_emax))
     211-                       mpfr_overflow (dest, rnd, sign) );
     212-  /* inexact = mpfr_set (dest, cache->x, rnd); */
     213-  if (MPFR_LIKELY(cache->inexact != 0))
     214+  MPFR_RNDRAW_GEN (inexact, dest,
     215+                   MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign,
     216+                   if (MPFR_UNLIKELY (cache->inexact == 0))
     217+                     {
     218+                       if ((sp[0] & ulp) == 0)
     219+                         {
     220+                           inexact = -sign;
     221+                           goto trunc_doit;
     222+                         }
     223+                       else
     224+                         goto addoneulp;
     225+                     }
     226+                   else if (cache->inexact < 0)
     227+                     goto addoneulp;
     228+                   else
     229+                     {
     230+                       inexact = -sign;
     231+                       goto trunc_doit;
     232+                     },
     233+                   if (MPFR_UNLIKELY (++MPFR_EXP (dest) > __gmpfr_emax))
     234+                     mpfr_overflow (dest, rnd, sign);
     235+                  );
     236+  if (MPFR_LIKELY (cache->inexact != 0))
     237     {
     238       switch (rnd)
     239         {
     240         case GMP_RNDZ:
     241         case GMP_RNDD:
     242-          if (MPFR_UNLIKELY(inexact == 0))
     243+          if (MPFR_UNLIKELY (inexact == 0))
     244             {
     245               inexact = cache->inexact;
     246               if (inexact > 0)
     247@@ -93,7 +106,7 @@
     248             }
     249           break;
     250         case GMP_RNDU:
     251-          if (MPFR_UNLIKELY(inexact == 0))
     252+          if (MPFR_UNLIKELY (inexact == 0))
     253             {
     254               inexact = cache->inexact;
     255               if (inexact < 0)
     256@@ -101,16 +114,7 @@
     257             }
     258           break;
     259         default: /* GMP_RNDN */
     260-          if (MPFR_UNLIKELY(inexact == MPFR_EVEN_INEX ||
     261-                            inexact == -MPFR_EVEN_INEX))
     262-            {
     263-              if (cache->inexact < 0)
     264-                mpfr_nextabove (dest);
     265-              else
     266-                mpfr_nextbelow (dest);
     267-              inexact = -inexact;
     268-            }
     269-          else if (MPFR_UNLIKELY(inexact == 0))
     270+          if (MPFR_UNLIKELY(inexact == 0))
     271             inexact = cache->inexact;
     272           break;
     273         }
     274diff -Naurd mpfr-2.2.0-p4/hypot.c mpfr-2.2.0-p5/hypot.c
     275--- hypot.c     2005-06-06 13:39:39.000000000 +0000
     276+++ hypot.c     2005-11-23 09:04:29.000000000 +0000
     277@@ -73,6 +73,7 @@
     278   Ey = MPFR_GET_EXP (y);
     279   diff_exp = (mp_exp_unsigned_t) Ex - Ey;
     280 
     281+  Nx = MPFR_PREC (x);   /* Precision of input variable */
     282   Nz = MPFR_PREC (z);   /* Precision of output variable */
     283 
     284   /* we have x < 2^Ex thus x^2 < 2^(2*Ex),
     285@@ -87,10 +88,10 @@
     286      if 2*diff_exp > Nx (see above as if Nz = Nx), therefore on Nz bits.
     287      Hence the condition: 2*diff_exp > MAX(Nz,Nx).
     288   */
     289-  if (diff_exp > MAX (Nz, MPFR_PREC (x)) / 2)
     290+  if (diff_exp > MAX (Nz, Nx) / 2)
     291     /* result is |x| or |x|+ulp(|x|,Nz) */
     292     {
     293-      if (rnd_mode == GMP_RNDU)
     294+      if (MPFR_UNLIKELY (rnd_mode == GMP_RNDU))
     295         {
     296           /* If z > abs(x), then it was already rounded up; otherwise
     297              z = abs(x), and we need to add one ulp due to y. */
     298@@ -100,14 +101,27 @@
     299         }
     300       else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */
     301         {
     302-          inexact = mpfr_abs (z, x, rnd_mode);
     303-          return (inexact) ? inexact : -1;
     304+          if (MPFR_LIKELY (Nz >= Nx))
     305+            {
     306+              mpfr_abs (z, x, rnd_mode);  /* exact */
     307+              return -1;
     308+            }
     309+          else
     310+            {
     311+              MPFR_SET_EXP (z, Ex);
     312+              MPFR_SET_SIGN (z, 1);
     313+              MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), Nx, rnd_mode, 1,
     314+                               goto addoneulp,
     315+                  if (MPFR_UNLIKELY (++MPFR_EXP (z) > __gmpfr_emax))
     316+                    return mpfr_overflow (z, rnd_mode, 1);
     317+                              );
     318+              return inexact ? inexact : -1;
     319+            }
     320         }
     321     }
     322 
     323   /* General case */
     324 
     325-  Nx = MPFR_PREC(x);   /* Precision of input variable */
     326   Ny = MPFR_PREC(y);   /* Precision of input variable */
     327 
     328   /* compute the working precision -- see algorithms.ps */
     329diff -Naurd mpfr-2.2.0-p4/mpfr-impl.h mpfr-2.2.0-p5/mpfr-impl.h
     330--- mpfr-impl.h 2005-09-11 22:13:24.000000000 +0000
     331+++ mpfr-impl.h 2005-11-23 09:04:29.000000000 +0000
     332@@ -184,6 +184,14 @@
     333 # define MPFR_THREAD_ATTR
     334 #endif
     335 
     336+/* Cache struct */
     337+struct __gmpfr_cache_s {
     338+  mpfr_t x;
     339+  int inexact;
     340+  int (*func)(mpfr_ptr, mpfr_rnd_t);
     341+};
     342+typedef struct __gmpfr_cache_s mpfr_cache_t[1];
     343+
     344 #if defined (__cplusplus)
     345 extern "C" {
     346 #endif
     347@@ -924,113 +932,19 @@
     348  ******************************************************/
     349 
     350 /*
     351- * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd`
     352- * assuming dest's sign is `sign`.
     353- * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case)
     354+ * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than
     355+ * once in a function (otherwise these labels would not be unique).
     356  */
     357-#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\
     358-  do {                                                                      \
     359-    mp_size_t dests, srcs;                                                  \
     360-    mp_limb_t *destp;                                                       \
     361-    mp_prec_t destprec, srcprec;                                            \
     362-                                                                            \
     363-    /* Check Trivial Case when Dest Mantissa has more bits than source */   \
     364-    srcprec = sprec;                                                        \
     365-    destprec = MPFR_PREC (dest);                                            \
     366-    destp = MPFR_MANT (dest);                                               \
     367-    if (MPFR_UNLIKELY (destprec >= srcprec))                                \
     368-      {                                                                     \
     369-        srcs  = (srcprec  + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
     370-        dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - srcs;    \
     371-        MPN_COPY (destp + dests, srcp, srcs);                               \
     372-        MPN_ZERO (destp, dests);                                            \
     373-        inexact = 0;                                                        \
     374-      }                                                                     \
     375-    else                                                                    \
     376-      {                                                                     \
     377-        /* Non trivial case: rounding needed */                             \
     378-        mp_prec_t sh;                                                       \
     379-        mp_limb_t *sp;                                                      \
     380-        mp_limb_t rb, sb, ulp;                                              \
     381-                                                                            \
     382-        /* Compute Position and shift */                                    \
     383-        srcs  = (srcprec  + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
     384-        dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
     385-        MPFR_UNSIGNED_MINUS_MODULO (sh, destprec);                          \
     386-        sp = srcp + srcs - dests;                                           \
     387-                                                                            \
     388-        /* General case when prec % BITS_PER_MP_LIMB != 0 */                \
     389-        if (MPFR_LIKELY (sh != 0))                                          \
     390-          {                                                                 \
     391-            mp_limb_t mask;                                                 \
     392-            /* Compute Rounding Bit and Sticky Bit */                       \
     393-            mask = MPFR_LIMB_ONE << (sh-1);                                 \
     394-            rb = sp[0] & mask;                                              \
     395-            sb = sp[0] & (mask-1);                                          \
     396-            if (MPFR_UNLIKELY (sb == 0))                                    \
     397-              { /* TODO: Improve it */                                      \
     398-                mp_limb_t *tmp;                                             \
     399-                mp_size_t n;                                                \
     400-                for (tmp = sp, n = srcs - dests ; n != 0 && sb == 0 ; n--)  \
     401-                  sb = *--tmp;                                              \
     402-              }                                                             \
     403-            ulp = 2*mask;                                                   \
     404-          }                                                                 \
     405-        else /* sh == 0 */                                                  \
     406-          {                                                                 \
     407-            MPFR_ASSERTD (dests < srcs);                                    \
     408-            /* Compute Rounding Bit and Sticky Bit */                       \
     409-            rb = sp[-1] & MPFR_LIMB_HIGHBIT;                                \
     410-            sb = sp[-1] & (MPFR_LIMB_HIGHBIT-1);                            \
     411-            if (MPFR_UNLIKELY (sb == 0))                                    \
     412-              {                                                             \
     413-                mp_limb_t *tmp;                                             \
     414-                mp_size_t n;                                                \
     415-                for (tmp = sp-1, n = srcs - dests-1 ; n!=0 && sb==0 ; n--)  \
     416-                  sb = *--tmp;                                              \
     417-              }                                                             \
     418-            ulp = MPFR_LIMB_ONE;                                            \
     419-          }                                                                 \
     420-        /* Rounding */                                                      \
     421-        if (MPFR_LIKELY (rnd == GMP_RNDN))                                  \
     422-          {                                                                 \
     423-            if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (sp[0] & ulp) == 0))   \
     424-              {                                                             \
     425-              trunc:                                                        \
     426-                inexact = MPFR_LIKELY ((sb | rb) != 0) ? -sign : 0;         \
     427-                MPN_COPY (destp, sp, dests);                                \
     428-                destp[0] &= ~(ulp-1);                                       \
     429-              }                                                             \
     430-            else                                                            \
     431-              {                                                             \
     432-              addoneulp:                                                    \
     433-                if (MPFR_UNLIKELY (mpn_add_1 (destp, sp, dests, ulp)))      \
     434-                  {                                                         \
     435-                    destp[dests-1] = MPFR_LIMB_HIGHBIT;                     \
     436-                    OVERFLOW_HANDLER;                                       \
     437-                  }                                                         \
     438-                destp[0] &= ~(ulp-1);                                       \
     439-                inexact = sign;                                             \
     440-              }                                                             \
     441-          }                                                                 \
     442-        else                                                                \
     443-          { /* Not Rounding to Nearest */                                   \
     444-            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\
     445-                || MPFR_UNLIKELY ((sb | rb) == 0))                          \
     446-              goto trunc;                                                   \
     447-             else                                                           \
     448-              goto addoneulp;                                               \
     449-          }                                                                 \
     450-      }                                                                     \
     451-  } while (0)
     452 
     453 /*
     454- * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd`
     455- * assuming dest's sign is `sign`.
     456- * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case)
     457- * Return MPFR_EVEN_INEX in case of EVEN rounding
     458+ * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
     459+ * assuming dest's sign is sign.
     460+ * In rounding to nearest mode, execute MIDDLE_HANDLER when the value
     461+ * is the middle of two consecutive numbers in dest precision.
     462+ * Execute OVERFLOW_HANDLER in case of overflow when rounding.
     463  */
     464-#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\
     465+#define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign,              \
     466+                        MIDDLE_HANDLER, OVERFLOW_HANDLER)                   \
     467   do {                                                                      \
     468     mp_size_t dests, srcs;                                                  \
     469     mp_limb_t *destp;                                                       \
     470@@ -1105,19 +1019,8 @@
     471                 destp[0] &= ~(ulp-1);                                       \
     472               }                                                             \
     473             else if (MPFR_UNLIKELY (sb == 0))                               \
     474-              {                                                             \
     475-                /* EVEN rounding */                                         \
     476-                if ((sp[0] & ulp) == 0)                                     \
     477-                 {                                                          \
     478-                  MPFR_ASSERTD (rb != 0);                                   \
     479-                  inexact = -MPFR_EVEN_INEX*sign;                           \
     480-                  goto trunc_doit;                                          \
     481-                 }                                                          \
     482-                else                                                        \
     483-                 {                                                          \
     484-                  inexact = MPFR_EVEN_INEX*sign;                            \
     485-                  goto addoneulp_doit;                                      \
     486-                 }                                                          \
     487+              { /* Middle of two consecutive representable numbers */       \
     488+                MIDDLE_HANDLER;                                             \
     489               }                                                             \
     490             else                                                            \
     491               {                                                             \
     492@@ -1133,16 +1036,58 @@
     493               }                                                             \
     494           }                                                                 \
     495         else                                                                \
     496-          { /* Not Rounding to Nearest */                                   \
     497-            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\
     498-                || MPFR_UNLIKELY ((sb | rb) == 0))                          \
     499+          { /* Directed rounding mode */                                    \
     500+            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd,                        \
     501+                                                MPFR_IS_NEG_SIGN (sign))))  \
     502               goto trunc;                                                   \
     503+             else if (MPFR_UNLIKELY ((sb | rb) == 0))                       \
     504+               {                                                            \
     505+                 inexact = 0;                                               \
     506+                 goto trunc_doit;                                           \
     507+               }                                                            \
     508              else                                                           \
     509               goto addoneulp;                                               \
     510           }                                                                 \
     511       }                                                                     \
     512   } while (0)
     513 
     514+/*
     515+ * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
     516+ * assuming dest's sign is sign.
     517+ * Execute OVERFLOW_HANDLER in case of overflow when rounding.
     518+ */
     519+#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \
     520+   MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,                   \
     521+     if ((sp[0] & ulp) == 0)                                                 \
     522+       {                                                                     \
     523+         inexact = -sign;                                                    \
     524+         goto trunc_doit;                                                    \
     525+       }                                                                     \
     526+     else                                                                    \
     527+       goto addoneulp;                                                       \
     528+     , OVERFLOW_HANDLER)
     529+
     530+/*
     531+ * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
     532+ * assuming dest's sign is sign.
     533+ * Execute OVERFLOW_HANDLER in case of overflow when rounding.
     534+ * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding.
     535+ */
     536+#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \
     537+                         OVERFLOW_HANDLER)                      \
     538+   MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,      \
     539+     if ((sp[0] & ulp) == 0)                                    \
     540+       {                                                        \
     541+         inexact = -MPFR_EVEN_INEX * sign;                      \
     542+         goto trunc_doit;                                       \
     543+       }                                                        \
     544+     else                                                       \
     545+       {                                                        \
     546+         inexact = MPFR_EVEN_INEX * sign;                       \
     547+         goto addoneulp_doit;                                   \
     548+       }                                                        \
     549+     , OVERFLOW_HANDLER)
     550+
     551 /* Return TRUE if b is non singular and we can round it to precision 'prec'
     552    with rounding mode 'rnd', and with error at most 'error' */
     553 #define MPFR_CAN_ROUND(b,err,prec,rnd)                                       \
     554@@ -1498,6 +1443,13 @@
     555 __MPFR_DECLSPEC int mpfr_const_log2_internal _MPFR_PROTO((mpfr_ptr,mp_rnd_t));
     556 __MPFR_DECLSPEC int mpfr_const_euler_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
     557 __MPFR_DECLSPEC int mpfr_const_catalan_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
     558+
     559+__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t,
     560+                                           int(*)(mpfr_ptr,mpfr_rnd_t)));
     561+__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t));
     562+__MPFR_DECLSPEC int  mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t,
     563+                                              mpfr_rnd_t));
     564+
     565 __MPFR_DECLSPEC void mpfr_mulhigh_n _MPFR_PROTO ((mp_ptr, mp_srcptr,
     566                                                   mp_srcptr, mp_size_t));
     567 
     568diff -Naurd mpfr-2.2.0-p4/mpfr.h mpfr-2.2.0-p5/mpfr.h
     569--- mpfr.h      2005-09-29 11:36:36.000000000 +0000
     570+++ mpfr.h      2005-11-23 09:04:29.000000000 +0000
     571@@ -118,14 +118,6 @@
     572 /* For those who needs a direct access and fast access to the sign field */
     573 #define MPFR_SIGN(x) ((x)->_mpfr_sign)
     574 
     575-/* Cache struct */
     576-struct __gmpfr_cache_s {
     577-  mpfr_t x;
     578-  int inexact;
     579-  int (*func)(mpfr_ptr, mpfr_rnd_t);
     580-};
     581-typedef struct __gmpfr_cache_s mpfr_cache_t[1];
     582-
     583 /* Stack interface */
     584 typedef enum {
     585   MPFR_NAN_KIND = 0,
     586@@ -544,11 +536,6 @@
     587 __MPFR_DECLSPEC int mpfr_sum _MPFR_PROTO ((mpfr_ptr, mpfr_ptr *__gmp_const,
     588                                            unsigned long, mpfr_rnd_t));
     589 
     590-__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t,
     591-                                           int(*)(mpfr_ptr,mpfr_rnd_t)));
     592-__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t));
     593-__MPFR_DECLSPEC int  mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t,
     594-                                              mpfr_rnd_t));
     595 __MPFR_DECLSPEC void mpfr_free_cache _MPFR_PROTO ((void));
     596 
     597 __MPFR_DECLSPEC int  mpfr_subnormalize _MPFR_PROTO ((mpfr_ptr, int,
     598diff -Naurd mpfr-2.2.0-p4/round_near_x.c mpfr-2.2.0-p5/round_near_x.c
     599--- round_near_x.c      2005-08-18 16:35:12.000000000 +0000
     600+++ round_near_x.c      2005-11-23 09:04:29.000000000 +0000
     601@@ -21,13 +21,13 @@
     602 
     603 #include "mpfr-impl.h"
     604 
     605-/* Uses MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
     606+/* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
     607 
     608 /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir,
     609                           mp_rnd_t rnd)
     610 
     611    Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(x)-error)
     612-   If x is small enought, y ~= x. This function checks and does this.
     613+   If x is small enough, y ~= x. This function checks and does this.
     614 
     615    It assumes that f(x) is not representable exactly as a FP number.
     616    x must not be a singular value (NAN, INF or ZERO).
     617@@ -42,7 +42,7 @@
     618    Otherwise it returns the ternary flag (It can't return an exact value).
     619 */
     620 
     621-/* What "small enought" means?
     622+/* What "small enough" means?
     623 
     624    We work with the positive values.
     625    Assuming err > Prec (y)+1
     626@@ -50,7 +50,7 @@
     627    i = [ y = o(x)]   // i = inexact flag
     628    If i == 0
     629        Setting x in y is exact. We have:
     630-       y = [XXXXXXXXX[...]]0[...] + error where [..] are optionnal zeros
     631+       y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros
     632       if dirError = ToInf,
     633         x < f(x) < x + 2^(EXP(x)-err)
     634         since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have:
     635@@ -172,9 +172,17 @@
     636   sign = MPFR_SIGN (x);
     637   MPFR_SET_EXP (y, MPFR_GET_EXP (x));
     638   MPFR_SET_SIGN (y, sign);
     639-  MPFR_RNDRAW_EVEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign,
     640-                    if (MPFR_UNLIKELY ( ++MPFR_EXP (y) > __gmpfr_emax))
     641-                    mpfr_overflow (y, rnd, sign) );
     642+  MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign,
     643+                   if (dir == 0)
     644+                     {
     645+                       inexact = -sign;
     646+                       goto trunc_doit;
     647+                     }
     648+                   else
     649+                     goto addoneulp;
     650+                   , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax))
     651+                       mpfr_overflow (y, rnd, sign)
     652+                  );
     653 
     654   /* Fix it in some cases */
     655   MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y));
     656@@ -212,15 +220,6 @@
     657             }
     658         }
     659     }
     660-  /* The even rule has been used. But due to error term, we should never
     661-     use this rule. That's why we have to fix some wrong rounding */
     662-  else if (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX)
     663-    {
     664-      if (inexact*sign > 0 && dir == 0)
     665-        goto nexttozero;
     666-      else if (inexact*sign < 0 && dir == 1)
     667-        goto nexttoinf;
     668-    }
     669 
     670   MPFR_RET (inexact);
     671 }
     672diff -Naurd mpfr-2.2.0-p4/tests/thypot.c mpfr-2.2.0-p5/tests/thypot.c
     673--- tests/thypot.c      2005-06-02 16:12:04.000000000 +0000
     674+++ tests/thypot.c      2005-11-23 09:04:29.000000000 +0000
     675@@ -131,7 +131,26 @@
     676   inexact = mpfr_hypot (z, x, y, GMP_RNDN);
     677   if (inexact >= 0 || mpfr_cmp (x, z))
     678     {
     679-      printf ("Error in test_large_small\n");
     680+      printf ("Error 1 in test_large_small\n");
     681+      exit (1);
     682+    }
     683+
     684+  mpfr_mul_ui (x, x, 5, GMP_RNDN);
     685+  inexact = mpfr_hypot (z, x, y, GMP_RNDN);
     686+  if (mpfr_cmp (x, z) >= 0)
     687+    {
     688+      printf ("Error 2 in test_large_small\n");
     689+      printf ("x = ");
     690+      mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
     691+      printf ("\n");
     692+      printf ("y = ");
     693+      mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
     694+      printf ("\n");
     695+      printf ("z = ");
     696+      mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
     697+      printf (" (in precision 2) instead of\n    ");
     698+      mpfr_out_str (stdout, 2, 2, x, GMP_RNDU);
     699+      printf ("\n");
     700       exit (1);
     701     }
     702 
     703diff -Naurd mpfr-2.2.0-p5/div.c mpfr-2.2.0-p6/div.c
     704--- div.c       2005-08-18 17:03:05.000000000 +0000
     705+++ div.c       2005-11-24 21:39:31.000000000 +0000
     706@@ -298,17 +298,16 @@
     707           MPN_COPY(bp, vp, vsize);
     708         }
     709       sticky_v = sticky_v || mpn_cmpzero (vp, k);
     710+      k = 0;
     711     }
     712-  else /* vsize < qsize */
     713+  else /* vsize < qsize: small divisor case */
     714     {
     715+      bp = vp;
     716       k = qsize - vsize;
     717-      bp = (mp_ptr) MPFR_TMP_ALLOC (qsize * sizeof(mp_limb_t));
     718-      MPN_COPY(bp + k, vp, vsize);
     719-      MPN_ZERO(bp, k);
     720     }
     721 
     722   /* we now can perform the division */
     723-  qh = mpn_divrem (qp, 0, ap, qqsize, bp, qsize);
     724+  qh = mpn_divrem (qp, 0, ap + k, qqsize - k, bp, qsize - k);
     725   /* warning: qh may be 1 if u1 == v1, but u < v */
     726 #ifdef DEBUG
     727   printf ("q="); mpn_print (qp, qsize);
     728diff -Naurd mpfr-2.2.0-p6/sin.c mpfr-2.2.0-p7/sin.c
     729--- sin.c       2005-08-18 17:03:11.000000000 +0000
     730+++ sin.c       2005-12-24 15:17:54.000000000 +0000
     731@@ -162,10 +162,12 @@
     732         {
     733           /* the absolute error on c is at most 2^(3-m-EXP(c)) */
     734           e = 2 * MPFR_GET_EXP (c) + m - 3;
     735-          if (mpfr_can_round (c, e, GMP_RNDZ, GMP_RNDZ,
     736+          if (mpfr_can_round (c, e, GMP_RNDN, GMP_RNDZ,
     737                               precy + (rnd_mode == GMP_RNDN)))
     738-            /* WARNING: need one more bit for rounding to nearest,
     739-               to be able to get the inexact flag correct */
     740+            /* WARNING: even if we know c <= sin(x), don't give GMP_RNDZ
     741+               as 3rd argument to mpfr_can_round, since if c is exactly
     742+               representable to the target precision (inexact = 0 below),
     743+               we would have to add one ulp when rounding away from 0. */
     744             break;
     745 
     746           /* check for huge cancellation (Near 0) */
     747@@ -183,14 +185,8 @@
     748   MPFR_ZIV_FREE (loop);
     749 
     750   inexact = mpfr_set (y, c, rnd_mode);
     751-
     752-  /* sin(x) is exact only for x = 0, which was treated apart above;
     753-     nevertheless, we can have inexact = 0 here if the approximation c
     754-     is exactly representable with PREC(y) bits. Since c is an approximation
     755-     towards zero, in that case the inexact flag should have the opposite sign
     756-     as y. */
     757-  if (MPFR_UNLIKELY (inexact == 0))
     758-    inexact = -MPFR_INT_SIGN (y);
     759+  /* inexact cannot be 0, since this would mean that c was representable
     760+     within the target precision, but in that case mpfr_can_round will fail */
     761 
     762   mpfr_clear (c);
     763 
Note: See TracChangeset for help on using the changeset viewer.