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

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