X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/rand/lcrand.c diff --git a/rand/lcrand.c b/rand/lcrand.c new file mode 100644 index 0000000..6944a71 --- /dev/null +++ b/rand/lcrand.c @@ -0,0 +1,304 @@ +/* -*-c-*- + * + * Simple linear congruential generator + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include +#include + +#include +#include + +#include "grand.h" +#include "lcrand.h" + +/*----- Magic numbers -----------------------------------------------------*/ + +/* --- The generator parameters --- */ + +#define P LCRAND_P /* Modulus */ +#define A LCRAND_A /* Multiplier (primitive mod @p@) */ +#define C LCRAND_C /* Additive constant */ + +/* --- Precomputed values for modular reduction --- */ + +#define D 5 /* %$p = 2^{32} - d$% */ + +/* --- Other useful bits --- */ + +#define P256 4294967040u /* Highest multiple of 256 < %$p$% */ + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @lcrand@ --- * + * + * Arguments: @uint32 x@ = seed value + * + * Returns: New state of the generator. + * + * Use: Steps the generator. Returns %$ax + c \bmod p$%. + */ + +uint32 lcrand(uint32 x) +{ + uint32 a[2], xx[2]; + uint32 yy[2]; + + /* --- Unpack things into the arrays --- */ + + a[0] = U16(A); a[1] = U16(A >> 16); + xx[0] = U16(x); xx[1] = U16(x >> 16); + + /* --- Multiply everything together --- * + * + * This is plain old long multiplication, although it looks a bit strange. + * I set up the top and bottom partial products directly where they're + * supposed to be. The cross terms I add together, with the low 16 bits in + * @q@ and the high 32 bits in @p@. These I then add into the product. + */ + + { + uint32 p, q; + + yy[0] = a[0] * xx[0]; + yy[1] = a[1] * xx[1]; + + p = a[0] * xx[1]; + q = p + a[1] * xx[0]; + p = ((q < p) << 16) + (q >> 16); + q = U16(q) << 16; + + q += yy[0]; + if (q < yy[0]) + p++; + else + p += (q >> 16) >> 16; + yy[0] = q; + + yy[1] += p; + } + + /* --- Now reduce mod p --- * + * + * I'm using shifts and adds to do the multiply step here. This needs to + * be changed if @D@ ever becomes something other than 5. + */ + +#if D != 5 +# error "Change shift sequence!" +#endif + + { + uint32 q; + + q = yy[1]; + x = yy[0]; + + while (q) { + uint32 y, z; + y = q >> 30; + z = q << 2; + z += q; + if (z < q) + y++; + else + y += (q >> 16) >> 16; + q = y; + x += z; + if (x < z || x > P) + x -= P; + } + } + + /* --- Now add on the constant --- */ + + x += C; + if (x < C || x >= P) + x -= P; + + /* --- Done --- */ + + return (x); +} + +/* --- @lcrand_range@ --- * + * + * Arguments: @uint32 *x@ = pointer to seed value (updated) + * @uint32 m@ = limit allowable + * + * Returns: A uniformly distributed pseudorandom integer in the interval + * %$[0, m)$%. + */ + +uint32 lcrand_range(uint32 *x, uint32 m) +{ + uint32 xx = *x; + uint32 r = P - P % m; + do xx = lcrand(xx); while (xx >= r); + *x = xx; + return (xx % m); +} + +/*----- Generic interface -------------------------------------------------*/ + +typedef struct gctx { + grand r; + uint32 x; +} gctx; + +static void gdestroy(grand *r) +{ + gctx *g = (gctx *)r; + DESTROY(g); +} + +static int gmisc(grand *r, unsigned op, ...) +{ + gctx *g = (gctx *)r; + va_list ap; + int rc = 0; + va_start(ap, op); + + switch (op) { + case GRAND_CHECK: + switch (va_arg(ap, unsigned)) { + case GRAND_CHECK: + case GRAND_SEEDINT: + case GRAND_SEEDUINT32: + case GRAND_SEEDRAND: + rc = 1; + break; + default: + rc = 0; + break; + } + break; + case GRAND_SEEDINT: + g->x = va_arg(ap, unsigned); + break; + case GRAND_SEEDUINT32: + g->x = va_arg(ap, uint32); + break; + case GRAND_SEEDRAND: { + grand *rr = va_arg(ap, grand *); + uint32 x; + do x = rr->ops->word(rr); while (x >= P || x == LCRAND_FIXEDPT); + g->x = x; + } break; + default: + GRAND_BADOP; + break; + } + + va_end(ap); + return (rc); +} + +static uint32 graw(grand *r) +{ + gctx *g = (gctx *)r; + g->x = lcrand(g->x); + return (g->x); +} + +static octet gbyte(grand *r) +{ + gctx *g = (gctx *)r; + uint32 x = g->x; + do x = lcrand(x); while (x >= P256); + g->x = x; + return (x / (P256 / 256)); +} + +static uint32 grange(grand *r, uint32 l) +{ + gctx *g = (gctx *)r; + return (lcrand_range(&g->x, l)); +} + +static const grand_ops gops = { + "lcrand", + LCRAND_P, 0, + gmisc, gdestroy, + graw, gbyte, grand_word, grange, grand_fill +}; + +/* --- @lcrand_create@ --- * + * + * Arguments: @uint32 x@ = initial seed + * + * Returns: Pointer to a generic generator. + * + * Use: Constructs a generic generator interface over a linear + * congruential generator. + */ + +grand *lcrand_create(uint32 x) +{ + gctx *g = CREATE(gctx); + g->r.ops = &gops; + g->x = x; + return (&g->r); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int verify(dstr *v) +{ + uint32 x = *(uint32 *)v[0].buf; + uint32 y = *(uint32 *)v[1].buf; + uint32 z = lcrand(x); + int ok = 1; + if (y != z) { + fprintf(stderr, + "\n*** lcrand failed. lcrand(%lu) = %lu, expected %lu\n", + (unsigned long)x, (unsigned long)z, (unsigned long)y); + ok = 0; + } + return (ok); +} + +static test_chunk tests[] = { + { "lcrand", verify, { &type_uint32, &type_uint32, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + test_run(argc, argv, tests, SRCDIR"/t/lcrand"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/