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

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