Thu Mar 14 06:01:07 1996 Roland McGrath <roland@charlie-brown.gnu.ai.mit.edu>
[kopensolaris-gnu/glibc.git] / stdlib / erand48_r.c
1 /* Copyright (C) 1995 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3 Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, August 1995.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Library General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13 Library General Public License for more details.
14
15 You should have received a copy of the GNU Library General Public
16 License along with the GNU C Library; see the file COPYING.LIB.  If
17 not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
18 Boston, MA 02111-1307, USA.  */
19
20 #include <stdlib.h>
21 #include "gmp.h"
22 #include "gmp-mparam.h"
23 #include <float.h>
24
25
26 /* Function to construct a floating point number from an MP integer
27    containing the fraction bits, a base 2 exponent, and a sign flag.  */
28 extern double __mpn_construct_double (mp_srcptr mpn, int exponent, int neg);
29
30 int
31 erand48_r (xsubi, buffer, result)
32      unsigned short int xsubi[3];
33      struct drand48_data *buffer;
34      double *result;
35 {
36   mp_limb mpn[(3 * sizeof (unsigned short int) + sizeof (mp_limb) - 1)
37               / sizeof (mp_limb)];
38
39   /* Compute next state.  */
40   if (__drand48_iterate (xsubi, buffer) < 0)
41     return -1;
42
43   /* Build a 48-bit mpn containing the 48 random bits.  */
44
45 #if BITS_PER_MP_LIMB == 64
46   mpn[0] = (xsubi[0] << 32) | (xsubi[1] << 16) | xsubi[2];
47 #elif BITS_PER_MP_LIMB == 32
48   mpn[0] = (xsubi[1] << 16) | xsubi[2];
49   mpn[1] = xsubi[0];
50 #else
51  #error "BITS_PER_MP_LIMB value not handled"
52 #endif
53
54   /* Shift them up so they are most significant bits of the fraction.  */
55   __mpn_lshift (mpn, mpn, sizeof mpn / sizeof mpn[0], DBL_MANT_DIG - 48);
56
57   /* Construct a positive double using those bits for the fractional part,
58      and a zero exponent so the resulting FP number is [0.0,1.0).  */
59   *result = __mpn_construct_double (mpn, 0, 0);
60
61   return 0;
62 }