[Mpfi-commits] r322 - trunk/mpfi/src

ptheveny at users.gforge.inria.fr ptheveny at users.gforge.inria.fr
Jeu 24 Juin 13:59:04 CEST 2010


Author: ptheveny
Date: 2010-06-24 13:59:04 +0200 (Thu, 24 Jun 2010)
New Revision: 322

Modified:
   trunk/mpfi/src/urandom.c
Log:
Improve mpfi_urandom: it previously always returned NaN for an infinite interval.

Modified: trunk/mpfi/src/urandom.c
===================================================================
--- trunk/mpfi/src/urandom.c	2010-06-24 11:52:49 UTC (rev 321)
+++ trunk/mpfi/src/urandom.c	2010-06-24 11:59:04 UTC (rev 322)
@@ -35,29 +35,57 @@
   mpfr_t diam, fact;
   prec = mpfr_get_prec (m);
 
-  if ( MPFI_NAN_P(y) ) {
+  if (MPFI_NAN_P(y)) {
     mpfr_set_nan (m);
+    return;
   }
 
+  if (mpfr_equal_p (&(y->left), &(y->right))) {
+    mpfr_set (m, &(y->left), MPFI_RNDD);
+    return;
+  }
+
   mpfr_init2 (diam, prec);
   mpfr_init2 (fact, prec);
 
   mpfi_diam_abs (diam, y);
-  mpfr_urandomb (fact, state);
+  mpfr_urandomb (fact, state); /* fact lies between 0 and 1 */
 
-  /* fact lies between 0 and 1, the picked point lies at a relative
-     distance "fact" of the left endpoint:  m = inf + (sup - inf)*fact  */
-  mpfr_mul (diam, diam, fact, MPFI_RNDD);
+  if (mpfr_cmp_ui (diam, 1) <= 0) {
+    /* the picked point lies at a relative distance "fact" of the left
+       endpoint: m = inf + (sup - inf) * fact  */
+    mpfr_mul (fact, fact, diam, MPFI_RNDD);
+    /* FIXME: because of possible cancelation, the random distribution is
+       not uniform among the floating-point numbers in y */
+    mpfr_add (m, &(y->left), fact, MPFI_RNDD);
+  }
+  else {
+    mp_exp_t e;
+    if (mpfr_cmp_abs (&(y->left), &(y->right)) < 0) {
+      e = mpfr_inf_p (&(y->right)) ? mpfr_get_emax ()
+        : mpfr_get_exp (&(y->right));
+    }
+    else {
+      e = mpfr_inf_p (&(y->left)) ? mpfr_get_emax ()
+        : mpfr_get_exp (&(y->left));
+    }
+    e += 1;
+    /* resize fact in [0, 2^e] where e = 1 + max{exp(left), exp(right)} */
+    mpfr_mul_2exp (fact, fact, e, MPFI_RNDD);
+    mpfr_set (m, &(y->left), MPFI_RNDD);
+    if (mpfr_inf_p (m)) {
+      mpfr_nextabove (m);
+    }
+    /* m may be outside y */
+    mpfr_add (m, m, fact, MPFI_RNDD);
+  }
+  mpfr_clear (fact);
+  mpfr_clear (diam);
 
-  mpfr_add (m, &(y->left), diam, MPFI_RNDD);
-  if (!mpfr_nan_p (m)) {  /* Ensure that m belongs to y */
-    if (mpfr_cmp (m, &(y->left)) < 0)
-      mpfr_set (m, &(y->left), MPFI_RNDU);
+  /* Ensure that m belongs to y (if the precision is sufficient) */
+  if (mpfr_cmp (m, &(y->left)) < 0)
+    mpfr_set (m, &(y->left), MPFI_RNDU);
 
-    if (mpfr_cmp (&(y->right), m) < 0)
-      mpfr_set (m, &(y->right), MPFI_RNDD);
-  }
-
-  mpfr_clear (diam);
-  mpfr_clear (fact);
+  if (mpfr_cmp (&(y->right), m) < 0)
+    mpfr_set (m, &(y->right), MPFI_RNDD);
 }




More information about the Mpfi-commits mailing list