gen/rand48.h

   1 /*
   2  * Copyright (c) 1993 Martin Birgmeier
   3  * All rights reserved.
   4  *
   5  * You may redistribute unmodified or modified versions of this source
   6  * code provided that the above copyright notice and this and the
   7  * following conditions are retained.
   8  *
   9  * This software is provided ``as is'', and comes with no warranties
  10  * of any kind. I shall in no event be liable for anything that happens
  11  * to anyone/anything when using this software.
  12  *
  13  * $FreeBSD: src/lib/libc/gen/rand48.h,v 1.2 2002/02/01 01:32:19 obrien Exp $
  14  */
  15
  16 #ifndef _RAND48_H_
  17 #define _RAND48_H_
  18
  19 #include <math.h>
  20 #include <stdlib.h>
  21
  22 #define RAND48_SEED_0   (0x330e)
  23 #define RAND48_SEED_1   (0xabcd)
  24 #define RAND48_SEED_2   (0x1234)
  25 #define RAND48_MULT_0   (0xe66d)
  26 #define RAND48_MULT_1   (0xdeec)
  27 #define RAND48_MULT_2   (0x0005)
  28 #define RAND48_ADD      (0x000b)
  29
  30 typedef unsigned long long uint48;
  31
  32 extern uint48 _rand48_seed;
  33 extern uint48 _rand48_mult;
  34 extern uint48 _rand48_add;
  35
  36 #define TOUINT48(x,y,z) \
  37         ((uint48)(x) + (((uint48)(y)) << 16) + (((uint48)(z)) << 32))
  38
  39 #define RAND48_SEED     TOUINT48(RAND48_SEED_0, RAND48_SEED_1, RAND48_SEED_2)
  40 #define RAND48_MULT     TOUINT48(RAND48_MULT_0, RAND48_MULT_1, RAND48_MULT_2)
  41
  42 #define LOADRAND48(l,x) \
  43         (l) = TOUINT48((x)[0], (x)[1], (x)[2])
  44
  45 #define STORERAND48(l,x) \
  46         (x)[0] = (unsigned short)(l); \
  47         (x)[1] = (unsigned short)((l) >> 16); \
  48         (x)[2] = (unsigned short)((l) >> 32)
  49
  50 #define _DORAND48(l) \
  51         (l) = (l) * _rand48_mult + _rand48_add
  52
  53 #define DORAND48(l,x) \
  54         LOADRAND48(l, x); \
  55         _DORAND48(l); \
  56         STORERAND48(l, x)
  57
  58 #include "fpmath.h"
  59
  60 /*
  61  * Optimization for speed: avoid int-to-double conversion.  Assume doubles
  62  * are IEEE-754 and insert the bits directly.  To normalize, the (1 << 52)
  63  * is the hidden bit, which the first set bit is shifted to.
  64  */
  65 #define ERAND48_BEGIN \
  66         union { \
  67                 union IEEEd2bits ieee; \
  68                 unsigned long long l; \
  69         } u; \
  70         int s
  71
  72 #define ERAND48_END(x) \
  73         u.l = ((x) & 0xffffffffffffULL); \
  74         if (u.l == 0) \
  75                 return 0.0; \
  76         u.l <<= 5; \
  77         for(s = 0; !(u.l & (1LL << 52)); s++, u.l <<= 1) {} \
  78         u.ieee.bits.exp = 1022 - s; \
  79         return u.ieee.d
  80
  81 #endif /* _RAND48_H_ */