chiark - git - mdw - storin/blob - lcrand.c

   1 /* -*-c-*-
   2  *
   3  * $Id: lcrand.c,v 1.2 2000/07/02 15:21:20 mdw Exp $
   4  *
   5  * Simple linear congruential generator
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  * (c) 2000 Mark Wooding
   9  */
  10
  11 /*----- Licensing notice --------------------------------------------------*
  12  *
  13  * Copyright (c) 2000 Mark Wooding
  14  * All rights reserved.
  15  *
  16  * Redistribution and use in source and binary forms, with or without
  17  * modification, are permitted provided that the following conditions are
  18  * met:
  19  *
  20  * 1. Redistributions of source code must retain the above copyright
  21  *    notice, this list of conditions and the following disclaimer.
  22  *
  23  * 2, Redistributions in binary form must reproduce the above copyright
  24  *    notice, this list of conditions and the following disclaimer in the
  25  *    documentation and/or other materials provided with the distribution.
  26  *
  27  * 3. The name of the authors may not be used to endorse or promote
  28  *    products derived from this software without specific prior written
  29  *    permission.
  30  *
  31  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED
  32  * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  33  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN
  34  * NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  35  * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  36  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  37  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  38  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  39  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  40  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  41  * POSSIBILITY OF SUCH DAMAGE.
  42  *
  43  * Instead of accepting the above terms, you may redistribute and/or modify
  44  * this software under the terms of either the GNU General Public License,
  45  * or the GNU Library General Public License, published by the Free
  46  * Software Foundation; either version 2 of the License, or (at your
  47  * option) any later version.
  48  */
  49
  50 /*----- Revision history --------------------------------------------------*
  51  *
  52  * $Log: lcrand.c,v $
  53  * Revision 1.2  2000/07/02 15:21:20  mdw
  54  * Fix licence text.
  55  *
  56  * Revision 1.1  2000/05/21 11:28:30  mdw
  57  * Initial check-in.
  58  *
  59  * --- Past lives (Catacomb) --- *
  60  *
  61  * Revision 1.2  1999/12/13 15:34:01  mdw
  62  * Add support for seeding from a generic pseudorandom source.
  63  *
  64  * Revision 1.1  1999/12/10 23:15:27  mdw
  65  * Noncryptographic random number generator.
  66  *
  67  */
  68
  69 /*----- Header files ------------------------------------------------------*/
  70
  71 #include <stdio.h>
  72 #include <stdlib.h>
  73 #include <string.h>
  74
  75 #include "bits.h"
  76 #include "lcrand.h"
  77
  78 /*----- Magic numbers -----------------------------------------------------*/
  79
  80 /* --- The generator parameters --- */
  81
  82 #define P LCRAND_P                      /* Modulus */
  83 #define A LCRAND_A                      /* Multiplier (primitive mod @p@) */
  84 #define C LCRAND_C                      /* Additive constant */
  85
  86 /* --- Precomputed values for modular reduction --- */
  87
  88 #define D 5                             /* %$p = 2^{32} - d$% */
  89
  90 /* --- Other useful bits --- */
  91
  92 #define P256 4294967040u                /* Highest multiple of 256 < %$p$% */
  93
  94 /*----- Main code ---------------------------------------------------------*/
  95
  96 /* --- @lcrand@ --- *
  97  *
  98  * Arguments:   @uint32 x@ = seed value
  99  *
 100  * Returns:     New state of the generator.
 101  *
 102  * Use:         Steps the generator.  Returns %$ax + c \bmod p$%.
 103  */
 104
 105 uint32 lcrand(uint32 x)
 106 {
 107   uint32 a[2], xx[2];
 108   uint32 yy[2];
 109
 110   /* --- Unpack things into the arrays --- */
 111
 112   a[0] = U16(A); a[1] = U16(A >> 16);
 113   xx[0] = U16(x); xx[1] = U16(x >> 16);
 114
 115   /* --- Multiply everything together --- *
 116    *
 117    * This is plain old long multiplication, although it looks a bit strange.
 118    * I set up the top and bottom partial products directly where they're
 119    * supposed to be.  The cross terms I add together, with the low 16 bits in
 120    * @q@ and the high 32 bits in @p@.  These I then add into the product.
 121    */
 122
 123   {
 124     uint32 p, q;
 125
 126     yy[0] = a[0] * xx[0];
 127     yy[1] = a[1] * xx[1];
 128
 129     p = a[0] * xx[1];
 130     q = p + a[1] * xx[0];
 131     p = ((q < p) << 16) + (q >> 16);
 132     q = U16(q) << 16;
 133
 134     q += yy[0];
 135     if (q < yy[0])
 136       p++;
 137     else
 138       p += (q >> 16) >> 16;
 139     yy[0] = q;
 140
 141     yy[1] += p;
 142   }
 143
 144   /* --- Now reduce mod p --- *
 145    *
 146    * I'm using shifts and adds to do the multiply step here.  This needs to
 147    * be changed if @D@ ever becomes something other than 5.
 148    */
 149
 150 #if D != 5
 151 #  error "Change shift sequence!"
 152 #endif
 153
 154   {
 155     uint32 q;
 156
 157     q = yy[1];
 158     x = yy[0];
 159
 160     while (q) {
 161       uint32 y, z;
 162       y = q >> 30;
 163       z = q << 2;
 164       z += q;
 165       if (z < q)
 166         y++;
 167       else
 168         y += (q >> 16) >> 16;
 169       q = y;
 170       x += z;
 171       if (x < z || x > P)
 172         x -= P;
 173     }
 174   }
 175
 176   /* --- Now add on the constant --- */
 177
 178   x += C;
 179   if (x < C || x >= P)
 180     x -= P;
 181
 182   /* --- Done --- */
 183
 184   return (x);
 185 }
 186
 187 /* --- @lcrand_range@ --- *
 188  *
 189  * Arguments:   @uint32 *x@ = pointer to seed value (updated)
 190  *              @uint32 m@ = limit allowable
 191  *
 192  * Returns:     A uniformly distributed pseudorandom integer in the interval
 193  *              %$[0, m)$%.
 194  */
 195
 196 uint32 lcrand_range(uint32 *x, uint32 m)
 197 {
 198   uint32 xx = *x;
 199   uint32 r = P - P % m;
 200   do xx = lcrand(xx); while (xx >= r);
 201   *x = xx;
 202   return (xx / (r / m));
 203 }
 204
 205 /*----- That's all, folks -------------------------------------------------*/