| 1 | /* Copyright (C) 1995, 1997 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 not, | 
|---|
| 17 | write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, | 
|---|
| 18 | Boston, MA 02111-1307, USA.  */ | 
|---|
| 19 |  | 
|---|
| 20 | #include <ieee754.h> | 
|---|
| 21 | #include <unixmac.h> | 
|---|
| 22 | #include <limits.h> | 
|---|
| 23 |  | 
|---|
| 24 |  | 
|---|
| 25 | int | 
|---|
| 26 | erand48_r (xsubi, buffer, result) | 
|---|
| 27 | unsigned short int xsubi[3]; | 
|---|
| 28 | struct drand48_data *buffer; | 
|---|
| 29 | double *result; | 
|---|
| 30 | { | 
|---|
| 31 | union ieee754_double temp; | 
|---|
| 32 |  | 
|---|
| 33 | /* Compute next state.  */ | 
|---|
| 34 | if (__drand48_iterate (xsubi, buffer) < 0) | 
|---|
| 35 | return -1; | 
|---|
| 36 |  | 
|---|
| 37 | /* Construct a positive double with the 48 random bits distributed over | 
|---|
| 38 | its fractional part so the resulting FP number is [0.0,1.0).  */ | 
|---|
| 39 |  | 
|---|
| 40 | #if USHRT_MAX == 65535 | 
|---|
| 41 | temp.ieee.negative = 0; | 
|---|
| 42 | temp.ieee.exponent = IEEE754_DOUBLE_BIAS; | 
|---|
| 43 | temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12); | 
|---|
| 44 | temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4); | 
|---|
| 45 | /* Please note the lower 4 bits of mantissa1 are always 0.  */ | 
|---|
| 46 | *result = temp.d - 1.0; | 
|---|
| 47 | #else | 
|---|
| 48 | # error Unsupported size of short int | 
|---|
| 49 | #endif | 
|---|
| 50 |  | 
|---|
| 51 | return 0; | 
|---|
| 52 | } | 
|---|