| version 1.1, 2000/09/09 14:12:19 |
version 1.1.1.2, 2003/08/25 16:06:07 |
|
|
| /* mpfr_random -- generate a random floating-point number |
/* mpfr_random -- generate a random floating-point number |
| |
|
| Copyright (C) 1999 PolKA project, Inria Lorraine and Loria |
Copyright 1999, 2001 Free Software Foundation, Inc. |
| |
|
| This file is part of the MPFR Library. |
This file is part of the MPFR Library. |
| |
|
| The MPFR Library is free software; you can redistribute it and/or modify |
The MPFR Library is free software; you can redistribute it and/or modify |
| it under the terms of the GNU Library General Public License as published by |
it under the terms of the GNU Lesser General Public License as published by |
| the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
| option) any later version. |
option) any later version. |
| |
|
| The MPFR Library is distributed in the hope that it will be useful, but |
The MPFR Library is distributed in the hope that it will be useful, but |
| WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
| License for more details. |
License for more details. |
| |
|
| You should have received a copy of the GNU Library General Public License |
You should have received a copy of the GNU Lesser General Public License |
| along with the MPFR Library; see the file COPYING.LIB. If not, write to |
along with the MPFR Library; see the file COPYING.LIB. If not, write to |
| the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
| MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
| Line 24 MA 02111-1307, USA. */ |
|
| Line 24 MA 02111-1307, USA. */ |
|
| #include "gmp-impl.h" |
#include "gmp-impl.h" |
| #include "longlong.h" |
#include "longlong.h" |
| #include "mpfr.h" |
#include "mpfr.h" |
| |
#include "mpfr-impl.h" |
| |
|
| /* Computes a random mpfr in [0, 1[ with precision PREC */ |
/* Computes a random mpfr in [0, 1[ with precision MPFR_PREC */ |
| |
|
| extern long random _PROTO((void)); |
|
| extern int srandom _PROTO((unsigned int)); |
|
| |
|
| /* extracted from GNU mpf */ |
|
| #if defined (__hpux) || defined (__alpha) |
|
| /* HPUX lacks random(). DEC OSF/1 1.2 random() returns a double. */ |
|
| #define random mrand48 |
|
| #define srandom srand48 |
|
| #endif |
|
| |
|
| void |
void |
| #if __STDC__ |
mpfr_random (mpfr_ptr x) |
| mpfr_random(mpfr_ptr x) |
|
| #else |
|
| mpfr_random(x) |
|
| mpfr_ptr x; |
|
| #endif |
|
| { |
{ |
| mp_limb_t *xp; unsigned long xn, i, cnt, prec=PREC(x); |
mp_limb_t *xp; unsigned long xn, cnt, prec = MPFR_PREC(x); |
| |
|
| xp = MANT(x); |
MPFR_CLEAR_FLAGS(x); |
| |
xp = MPFR_MANT(x); |
| xn = (prec-1)/BITS_PER_MP_LIMB + 1; |
xn = (prec-1)/BITS_PER_MP_LIMB + 1; |
| |
|
| for (i = 0; i < xn; i++) |
mpn_random (xp, xn); |
| { |
|
| /* random() c/sh/ould be replaced by a homemade random number generator. |
|
| Indeed, if on Linux random is a good RNG, this is definitely not |
|
| the case in most Un*xes. */ |
|
| xp[i] = random(); |
|
| } |
|
| |
|
| count_leading_zeros(cnt, xp[xn - 1]); |
|
| if (cnt) mpn_lshift(xp, xp, xn, cnt); |
|
| EXP(x) = -cnt; |
|
| |
|
| cnt = xn*BITS_PER_MP_LIMB - prec; |
if (xp[xn - 1] == 0) |
| /* cnt is the number of non significant bits in the low limb */ |
xp[xn - 1] = 1; |
| xp[0] &= ~((((mp_limb_t)1)<<cnt) - 1); |
/* since count_leading_zeros doesn't like zeroes */ |
| } |
|
| |
|
| void |
count_leading_zeros (cnt, xp[xn - 1]); |
| mpfr_srandom(unsigned long seed) |
if (cnt) |
| { |
mpn_lshift (xp, xp, xn, cnt); |
| srandom(seed); |
MPFR_EXP(x) = -cnt; |
| |
if (MPFR_SIGN(x) < 0) |
| |
MPFR_CHANGE_SIGN(x); |
| |
|
| |
cnt = xn * BITS_PER_MP_LIMB - prec; |
| |
/* cnt is the number of non significant bits in the low limb */ |
| |
xp[0] &= ~((MP_LIMB_T_ONE << cnt) - MP_LIMB_T_ONE); |
| } |
} |