git.distorted.org.uk Git - u/mdw/catacomb/blob - lcrand.c

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