1 /* erandom.c -- Random element in the interval, following an exponential distribution.
2 without any guarantee: to be checked
3
4 Copyright 2018
5 AriC project, Inria Rhone-Alpes, France
6
7 This file is part of the MPFI Library.
8
9 The MPFI Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 2.1 of the License, or (at your
12 option) any later version.
13
14 The MPFI Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the MPFI Library; see the file COPYING.LIB. If not, write to
21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
22 MA 02110-1301, USA. */
23
24 #include "mpfi-impl.h"
25
26 /* Picks randomly a point m in y */
27 void
mpfi_erandom(mpfr_ptr m,mpfi_srcptr y,gmp_randstate_t state)28 mpfi_erandom (mpfr_ptr m, mpfi_srcptr y, gmp_randstate_t state)
29 {
30 mpfr_prec_t prec, tmp_prec;
31 mpfr_t diam, fact;
32
33
34 if (MPFI_NAN_P(y)) {
35 mpfr_set_nan (m);
36 return;
37 }
38
39 if (mpfr_equal_p (&(y->left), &(y->right))) {
40 mpfr_set (m, &(y->left), MPFR_RNDN);
41 return;
42 }
43
44 prec = mpfr_get_prec (m);
45 tmp_prec = mpfi_get_prec(y);
46 if (tmp_prec > prec)
47 {
48 prec = tmp_prec;
49 }
50 mpfr_init2 (diam, prec);
51 mpfr_init2 (fact, prec);
52
53 mpfi_diam_abs (diam, y);
54 mpfr_erandom (fact, state, MPFR_RNDN);
55 mpfr_sub_d (fact, fact, 0.5, MPFR_RNDN);
56 if (mpfr_cmp_ui (fact, 0) < 0)
57 mpfr_set_ui(fact, 0, MPFR_RNDN);
58 else if (mpfr_cmp_ui(fact, 1) > 0)
59 mpfr_set_ui(fact, 1, MPFR_RNDN); /* now fact lies between 0 and 1 */
60
61 if (mpfr_cmp_ui (diam, 1) <= 0) {
62 /* the picked point lies at a relative distance "fact" of the left
63 endpoint: m = inf + (sup - inf) * fact */
64 mpfr_mul (fact, fact, diam, MPFR_RNDN);
65 /* FIXME: because of possible cancelation, the random distribution is
66 not uniform among the floating-point numbers in y */
67 mpfr_add (m, &(y->left), fact, MPFR_RNDN);
68 }
69 else {
70 mpfr_exp_t e;
71 if (mpfr_cmp_abs (&(y->left), &(y->right)) < 0) {
72 e = mpfr_inf_p (&(y->right)) ? mpfr_get_emax ()
73 : mpfr_get_exp (&(y->right));
74 }
75 else {
76 e = mpfr_inf_p (&(y->left)) ? mpfr_get_emax ()
77 : mpfr_get_exp (&(y->left));
78 }
79 e += 1;
80 /* resize fact in [0, 2^e] where e = 1 + max{exp(left), exp(right)} */
81 mpfr_mul_2exp (fact, fact, e, MPFR_RNDN);
82 mpfr_set (m, &(y->left), MPFR_RNDN);
83 if (mpfr_inf_p (m)) {
84 mpfr_nextabove (m);
85 }
86 /* m may be outside y */
87 mpfr_add (m, m, fact, MPFR_RNDN);
88 }
89 mpfr_clear (fact);
90 mpfr_clear (diam);
91
92 /* Ensure that m belongs to y (if the precision is sufficient) */
93 if (mpfr_cmp (m, &(y->left)) < 0)
94 mpfr_set (m, &(y->left), MPFR_RNDU);
95
96 if (mpfr_cmp (&(y->right), m) < 0)
97 mpfr_set (m, &(y->right), MPFR_RNDD);
98 }
99