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

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