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

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