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