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

   1 /* -*-c-*-
   2  *
   3  * $Id: lcrand.c,v 1.3 2000/06/17 11:29:03 mdw Exp $
   4  *
   5  * Simple linear congruential generator
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: lcrand.c,v $
  33  * Revision 1.3  2000/06/17 11:29:03  mdw
  34  * Add the flags word to the generic generator.
  35  *
  36  * Revision 1.2  1999/12/13 15:34:01  mdw
  37  * Add support for seeding from a generic pseudorandom source.
  38  *
  39  * Revision 1.1  1999/12/10 23:15:27  mdw
  40  * Noncryptographic random number generator.
  41  *
  42  */
  43
  44 /*----- Header files ------------------------------------------------------*/
  45
  46 #include <stdarg.h>
  47 #include <stdio.h>
  48 #include <stdlib.h>
  49 #include <string.h>
  50
  51 #include <mLib/bits.h>
  52 #include <mLib/sub.h>
  53
  54 #include "grand.h"
  55 #include "lcrand.h"
  56
  57 /*----- Magic numbers -----------------------------------------------------*/
  58
  59 /* --- The generator parameters --- */
  60
  61 #define P LCRAND_P                      /* Modulus */
  62 #define A LCRAND_A                      /* Multiplier (primitive mod @p@) */
  63 #define C LCRAND_C                      /* Additive constant */
  64
  65 /* --- Precomputed values for modular reduction --- */
  66
  67 #define D 5                             /* %$p = 2^{32} - d$% */
  68
  69 /* --- Other useful bits --- */
  70
  71 #define P256 4294967040u                /* Highest multiple of 256 < %$p$% */
  72
  73 /*----- Main code ---------------------------------------------------------*/
  74
  75 /* --- @lcrand@ --- *
  76  *
  77  * Arguments:   @uint32 x@ = seed value
  78  *
  79  * Returns:     New state of the generator.
  80  *
  81  * Use:         Steps the generator.  Returns %$ax + c \bmod p$%.
  82  */
  83
  84 uint32 lcrand(uint32 x)
  85 {
  86   uint32 a[2], xx[2];
  87   uint32 yy[2];
  88
  89   /* --- Unpack things into the arrays --- */
  90
  91   a[0] = U16(A); a[1] = U16(A >> 16);
  92   xx[0] = U16(x); xx[1] = U16(x >> 16);
  93
  94   /* --- Multiply everything together --- *
  95    *
  96    * This is plain old long multiplication, although it looks a bit strange.
  97    * I set up the top and bottom partial products directly where they're
  98    * supposed to be.  The cross terms I add together, with the low 16 bits in
  99    * @q@ and the high 32 bits in @p@.  These I then add into the product.
 100    */
 101
 102   {
 103     uint32 p, q;
 104
 105     yy[0] = a[0] * xx[0];
 106     yy[1] = a[1] * xx[1];
 107
 108     p = a[0] * xx[1];
 109     q = p + a[1] * xx[0];
 110     p = ((q < p) << 16) + (q >> 16);
 111     q = U16(q) << 16;
 112
 113     q += yy[0];
 114     if (q < yy[0])
 115       p++;
 116     else
 117       p += (q >> 16) >> 16;
 118     yy[0] = q;
 119
 120     yy[1] += p;
 121   }
 122
 123   /* --- Now reduce mod p --- *
 124    *
 125    * I'm using shifts and adds to do the multiply step here.  This needs to
 126    * be changed if @D@ ever becomes something other than 5.
 127    */
 128
 129 #if D != 5
 130 #  error "Change shift sequence!"
 131 #endif
 132
 133   {
 134     uint32 q;
 135
 136     q = yy[1];
 137     x = yy[0];
 138
 139     while (q) {
 140       uint32 y, z;
 141       y = q >> 30;
 142       z = q << 2;
 143       z += q;
 144       if (z < q)
 145         y++;
 146       else
 147         y += (q >> 16) >> 16;
 148       q = y;
 149       x += z;
 150       if (x < z || x > P)
 151         x -= P;
 152     }
 153   }
 154
 155   /* --- Now add on the constant --- */
 156
 157   x += C;
 158   if (x < C || x >= P)
 159     x -= P;
 160
 161   /* --- Done --- */
 162
 163   return (x);
 164 }
 165
 166 /* --- @lcrand_range@ --- *
 167  *
 168  * Arguments:   @uint32 *x@ = pointer to seed value (updated)
 169  *              @uint32 m@ = limit allowable
 170  *
 171  * Returns:     A uniformly distributed pseudorandom integer in the interval
 172  *              %$[0, m)$%.
 173  */
 174
 175 uint32 lcrand_range(uint32 *x, uint32 m)
 176 {
 177   uint32 xx = *x;
 178   uint32 r = P - P % m;
 179   do xx = lcrand(xx); while (xx >= r);
 180   *x = xx;
 181   return (xx / (r / m));
 182 }
 183
 184 /*----- Generic interface -------------------------------------------------*/
 185
 186 typedef struct gctx {
 187   grand r;
 188   uint32 x;
 189 } gctx;
 190
 191 static void gdestroy(grand *r)
 192 {
 193   gctx *g = (gctx *)r;
 194   DESTROY(g);
 195 }
 196
 197 static int gmisc(grand *r, unsigned op, ...)
 198 {
 199   gctx *g = (gctx *)r;
 200   va_list ap;
 201   int rc = 0;
 202   va_start(ap, op);
 203
 204   switch (op) {
 205     case GRAND_CHECK:
 206       switch (va_arg(ap, unsigned)) {
 207         case GRAND_CHECK:
 208         case GRAND_SEEDINT:
 209         case GRAND_SEEDUINT32:
 210         case GRAND_SEEDRAND:
 211           rc = 1;
 212           break;
 213         default:
 214           rc = 0;
 215           break;
 216       }
 217       break;
 218     case GRAND_SEEDINT:
 219       g->x = va_arg(ap, unsigned);
 220       break;
 221     case GRAND_SEEDUINT32:
 222       g->x = va_arg(ap, uint32);
 223       break;
 224     case GRAND_SEEDRAND: {
 225       grand *rr = va_arg(ap, grand *);
 226       uint32 x;
 227       do x = rr->ops->word(rr); while (x >= P || x == LCRAND_FIXEDPT);
 228       g->x = x;
 229     } break;
 230     default:
 231       GRAND_BADOP;
 232       break;
 233   }
 234
 235   va_end(ap);
 236   return (rc);
 237 }
 238
 239 static uint32 graw(grand *r)
 240 {
 241   gctx *g = (gctx *)r;
 242   g->x = lcrand(g->x);
 243   return (g->x);
 244 }
 245
 246 static octet gbyte(grand *r)
 247 {
 248   gctx *g = (gctx *)r;
 249   uint32 x = g->x;
 250   do x = lcrand(x); while (x >= P256);
 251   g->x = x;
 252   return (x / (P256 / 256));
 253 }
 254
 255 static uint32 grange(grand *r, uint32 l)
 256 {
 257   gctx *g = (gctx *)r;
 258   return (lcrand_range(&g->x, l));
 259 }
 260
 261 static const grand_ops gops = {
 262   "lcrand",
 263   LCRAND_P, 0,
 264   gmisc, gdestroy,
 265   graw, gbyte, grand_word, grange, grand_fill
 266 };
 267
 268 /* --- @lcrand_create@ --- *
 269  *
 270  * Arguments:   @uint32 x@ = initial seed
 271  *
 272  * Returns:     Pointer to a generic generator.
 273  *
 274  * Use:         Constructs a generic generator interface over a linear
 275  *              congruential generator.
 276  */
 277
 278 grand *lcrand_create(uint32 x)
 279 {
 280   gctx *g = CREATE(gctx);
 281   g->r.ops = &gops;
 282   g->x = x;
 283   return (&g->r);
 284 }
 285
 286 /*----- Test rig ----------------------------------------------------------*/
 287
 288 #ifdef TEST_RIG
 289
 290 #include <mLib/testrig.h>
 291
 292 static int verify(dstr *v)
 293 {
 294   uint32 x = *(uint32 *)v[0].buf;
 295   uint32 y = *(uint32 *)v[1].buf;
 296   uint32 z = lcrand(x);
 297   int ok = 1;
 298   if (y != z) {
 299     fprintf(stderr,
 300             "\n*** lcrand failed.  lcrand(%lu) = %lu, expected %lu\n",
 301             (unsigned long)x, (unsigned long)z, (unsigned long)y);
 302     ok = 0;
 303   }
 304   return (ok);
 305 }
 306
 307 static test_chunk tests[] = {
 308   { "lcrand", verify, { &type_uint32, &type_uint32, 0 } },
 309   { 0, 0, { 0 } }
 310 };
 311
 312 int main(int argc, char *argv[])
 313 {
 314   test_run(argc, argv, tests, SRCDIR"/tests/lcrand");
 315   return (0);
 316 }
 317
 318 #endif
 319
 320 /*----- That's all, folks -------------------------------------------------*/