summaryrefslogtreecommitdiff
path: root/libm/ldouble/ldrand.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/ldouble/ldrand.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/ldouble/ldrand.c')
-rw-r--r--libm/ldouble/ldrand.c175
1 files changed, 175 insertions, 0 deletions
diff --git a/libm/ldouble/ldrand.c b/libm/ldouble/ldrand.c
new file mode 100644
index 000000000..892b465df
--- /dev/null
+++ b/libm/ldouble/ldrand.c
@@ -0,0 +1,175 @@
+/* ldrand.c
+ *
+ * Pseudorandom number generator
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double y;
+ * int ldrand();
+ *
+ * ldrand( &y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Yields a random number 1.0 <= y < 2.0.
+ *
+ * The three-generator congruential algorithm by Brian
+ * Wichmann and David Hill (BYTE magazine, March, 1987,
+ * pp 127-8) is used.
+ *
+ * Versions invoked by the different arithmetic compile
+ * time options IBMPC, and MIEEE, produce the same sequences.
+ *
+ */
+
+
+
+#include <math.h>
+#ifdef ANSIPROT
+int ranwh ( void );
+#else
+int ranwh();
+#endif
+#ifdef UNK
+#undef UNK
+#if BIGENDIAN
+#define MIEEE
+#else
+#define IBMPC
+#endif
+#endif
+
+/* Three-generator random number algorithm
+ * of Brian Wichmann and David Hill
+ * BYTE magazine, March, 1987 pp 127-8
+ *
+ * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
+ */
+
+static int sx = 1;
+static int sy = 10000;
+static int sz = 3000;
+
+static union {
+ long double d;
+ unsigned short s[8];
+} unkans;
+
+/* This function implements the three
+ * congruential generators.
+ */
+
+int ranwh()
+{
+int r, s;
+
+/* sx = sx * 171 mod 30269 */
+r = sx/177;
+s = sx - 177 * r;
+sx = 171 * s - 2 * r;
+if( sx < 0 )
+ sx += 30269;
+
+
+/* sy = sy * 172 mod 30307 */
+r = sy/176;
+s = sy - 176 * r;
+sy = 172 * s - 35 * r;
+if( sy < 0 )
+ sy += 30307;
+
+/* sz = 170 * sz mod 30323 */
+r = sz/178;
+s = sz - 178 * r;
+sz = 170 * s - 63 * r;
+if( sz < 0 )
+ sz += 30323;
+/* The results are in static sx, sy, sz. */
+return 0;
+}
+
+/* ldrand.c
+ *
+ * Random double precision floating point number between 1 and 2.
+ *
+ * C callable:
+ * drand( &x );
+ */
+
+int ldrand( a )
+long double *a;
+{
+unsigned short r;
+
+/* This algorithm of Wichmann and Hill computes a floating point
+ * result:
+ */
+ranwh();
+unkans.d = sx/30269.0L + sy/30307.0L + sz/30323.0L;
+r = unkans.d;
+unkans.d -= r;
+unkans.d += 1.0L;
+
+if( sizeof(long double) == 16 )
+ {
+#ifdef MIEEE
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[7] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[6] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[5] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[4] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[3] = r;
+#endif
+#ifdef IBMPC
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[0] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[1] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[2] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[3] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[4] = r;
+#endif
+ }
+else
+ {
+#ifdef MIEEE
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[5] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[4] = r;
+#endif
+#ifdef IBMPC
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[0] = r;
+ ranwh();
+ r = sx * sy + sz;
+ unkans.s[1] = r;
+#endif
+ }
+*a = unkans.d;
+return 0;
+}