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

   1 /* -*-c-*-
   2  *
   3  * $Id: lcrand.c,v 1.1 2000/05/21 11:28:30 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 REGENTS 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.1  2000/05/21 11:28:30  mdw
  54  * Initial check-in.
  55  *
  56  * --- Past lives (Catacomb) --- *
  57  *
  58  * Revision 1.2  1999/12/13 15:34:01  mdw
  59  * Add support for seeding from a generic pseudorandom source.
  60  *
  61  * Revision 1.1  1999/12/10 23:15:27  mdw
  62  * Noncryptographic random number generator.
  63  *
  64  */
  65
  66 /*----- Header files ------------------------------------------------------*/
  67
  68 #include <stdio.h>
  69 #include <stdlib.h>
  70 #include <string.h>
  71
  72 #include "bits.h"
  73 #include "lcrand.h"
  74
  75 /*----- Magic numbers -----------------------------------------------------*/
  76
  77 /* --- The generator parameters --- */
  78
  79 #define P LCRAND_P                      /* Modulus */
  80 #define A LCRAND_A                      /* Multiplier (primitive mod @p@) */
  81 #define C LCRAND_C                      /* Additive constant */
  82
  83 /* --- Precomputed values for modular reduction --- */
  84
  85 #define D 5                             /* %$p = 2^{32} - d$% */
  86
  87 /* --- Other useful bits --- */
  88
  89 #define P256 4294967040u                /* Highest multiple of 256 < %$p$% */
  90
  91 /*----- Main code ---------------------------------------------------------*/
  92
  93 /* --- @lcrand@ --- *
  94  *
  95  * Arguments:   @uint32 x@ = seed value
  96  *
  97  * Returns:     New state of the generator.
  98  *
  99  * Use:         Steps the generator.  Returns %$ax + c \bmod p$%.
 100  */
 101
 102 uint32 lcrand(uint32 x)
 103 {
 104   uint32 a[2], xx[2];
 105   uint32 yy[2];
 106
 107   /* --- Unpack things into the arrays --- */
 108
 109   a[0] = U16(A); a[1] = U16(A >> 16);
 110   xx[0] = U16(x); xx[1] = U16(x >> 16);
 111
 112   /* --- Multiply everything together --- *
 113    *
 114    * This is plain old long multiplication, although it looks a bit strange.
 115    * I set up the top and bottom partial products directly where they're
 116    * supposed to be.  The cross terms I add together, with the low 16 bits in
 117    * @q@ and the high 32 bits in @p@.  These I then add into the product.
 118    */
 119
 120   {
 121     uint32 p, q;
 122
 123     yy[0] = a[0] * xx[0];
 124     yy[1] = a[1] * xx[1];
 125
 126     p = a[0] * xx[1];
 127     q = p + a[1] * xx[0];
 128     p = ((q < p) << 16) + (q >> 16);
 129     q = U16(q) << 16;
 130
 131     q += yy[0];
 132     if (q < yy[0])
 133       p++;
 134     else
 135       p += (q >> 16) >> 16;
 136     yy[0] = q;
 137
 138     yy[1] += p;
 139   }
 140
 141   /* --- Now reduce mod p --- *
 142    *
 143    * I'm using shifts and adds to do the multiply step here.  This needs to
 144    * be changed if @D@ ever becomes something other than 5.
 145    */
 146
 147 #if D != 5
 148 #  error "Change shift sequence!"
 149 #endif
 150
 151   {
 152     uint32 q;
 153
 154     q = yy[1];
 155     x = yy[0];
 156
 157     while (q) {
 158       uint32 y, z;
 159       y = q >> 30;
 160       z = q << 2;
 161       z += q;
 162       if (z < q)
 163         y++;
 164       else
 165         y += (q >> 16) >> 16;
 166       q = y;
 167       x += z;
 168       if (x < z || x > P)
 169         x -= P;
 170     }
 171   }
 172
 173   /* --- Now add on the constant --- */
 174
 175   x += C;
 176   if (x < C || x >= P)
 177     x -= P;
 178
 179   /* --- Done --- */
 180
 181   return (x);
 182 }
 183
 184 /* --- @lcrand_range@ --- *
 185  *
 186  * Arguments:   @uint32 *x@ = pointer to seed value (updated)
 187  *              @uint32 m@ = limit allowable
 188  *
 189  * Returns:     A uniformly distributed pseudorandom integer in the interval
 190  *              %$[0, m)$%.
 191  */
 192
 193 uint32 lcrand_range(uint32 *x, uint32 m)
 194 {
 195   uint32 xx = *x;
 196   uint32 r = P - P % m;
 197   do xx = lcrand(xx); while (xx >= r);
 198   *x = xx;
 199   return (xx / (r / m));
 200 }
 201
 202 /*----- That's all, folks -------------------------------------------------*/