diff options
Diffstat (limited to 'rand.c')
| -rw-r--r-- | rand.c | 199 |
1 files changed, 199 insertions, 0 deletions
@@ -0,0 +1,199 @@ +/* Copyright 2012 + * Kaz Kylheku <kaz@kylheku.com> + * Vancouver, Canada + * All rights reserved. + * + * BSD License: + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * 3. The name of the author may not be used to endorse or promote + * products derived from this software without specific prior + * written permission. + * + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <wctype.h> +#include <assert.h> +#include <limits.h> +#include <stdarg.h> +#include <dirent.h> +#include <setjmp.h> +#include <wchar.h> +#include <limits.h> +#include <time.h> +#include "config.h" +#include "lib.h" +#include "unwind.h" +#include "gc.h" +#include "arith.h" +#include "rand.h" + +#if SIZEOF_INT == 4 +typedef unsigned int rand32_t; +#elif SIZEOF_LONG == 4 +typedef unsigned long rand32_t; +#endif + +struct random_state { + rand32_t state[16]; + int cur; +}; + +val random_state; +val random_state_s; + +static struct cobj_ops random_state_ops = { + cobj_equal_op, + cobj_print_op, + cobj_destroy_free_op, + cobj_mark_op, + cobj_hash_op +}; + +static val make_state() +{ + struct random_state *r = (struct random_state *) chk_malloc(sizeof *r); + return cobj((mem_t *) r, random_state_s, &random_state_ops); +} + +val random_state_p(val obj) +{ + return typeof(obj) == random_state_s ? t : nil; +} + +val make_random_state(val seed) +{ + val rs = make_state(); + struct random_state *r = (struct random_state *) + cobj_handle(rs, random_state_s); + + r->cur = 0; + + if (bignump(seed)) { + int i, dig, bit; + mp_int *m = mp(seed); + + for (i = 0, dig = 0, bit = 0; i < 16 && dig < m->used; i++) { + r->state[i] = (m->dp[dig] >> bit) & 0xFFFFFFFFul; + bit += 32; + if (bit >= MP_DIGIT_BIT) + dig++, bit = 0; + } + + for (; i < 16; i++) { + r->state[i] = 0xAAAAAAAA; + } + } else if (fixnump(seed)) { + cnum s = c_num(seed) & NUM_MAX; + + memset(r->state, 0xAA, sizeof r->state); + r->state[0] = s & 0xFFFFFFFFul; +#if SIZEOF_PTR >= 8 + r->state[1] = (s >> 32) & 0xFFFFFFFFul; +#elif SIZEOF_PTR >= 16 + r->state[2] = (s >> 64) & 0xFFFFFFFFul; + r->state[3] = (s >> 96) & 0xFFFFFFFFul; +#endif + } else if (nullp(seed)) { + r->state[0] = (rand32_t) time(0); + r->state[1] = (rand32_t) clock(); + memset(r->state + 2, 0xAA, sizeof r->state - 2 * sizeof r->state[0]); + } else if (random_state_p(seed)) { + struct random_state *rseed = (struct random_state *) + cobj_handle(seed, random_state_s); + *r = *rseed; + } else { + uw_throwf(error_s, lit("make-random-state: seed ~s is not a number"), + seed, nao); + } + + return rs; +} + +static rand32_t rand32(struct random_state *r) +{ + #define RSTATE(r,i) ((r)->state[((r)->cur + i) % 16]) + rand32_t s0 = RSTATE(r, 0); + rand32_t s9 = RSTATE(r, 9); + rand32_t s13 = RSTATE(r, 13); + rand32_t s15 = RSTATE(r, 15); + + rand32_t r1 = s0 ^ (s0 << 16) ^ s13 ^ (s13 << 15); + rand32_t r2 = s9 ^ (s9 >> 11); + + rand32_t ns0 = RSTATE(r, 0) = r1 ^ r2; + rand32_t ns15 = s15 ^ (s15 << 2) ^ r1 ^ (r1 << 18) ^ r2 ^ (r2 << 28) ^ + ((ns0 ^ (ns0 << 5)) & 0xda442d24ul); + + RSTATE(r, 15) = ns15; + r->cur = (r->cur + 15) % 16; + return ns15; + #undef RSTATE +} + +val random_fixnum(val state) +{ + uses_or2; + struct random_state *r = (struct random_state *) + cobj_handle(or2(state, random_state), + random_state_s); + return num(rand32(r) & NUM_MAX); +} + +val random(val state, val modulus) +{ + uses_or2; + struct random_state *r = (struct random_state *) + cobj_handle(or2(state, random_state), + random_state_s); + if (bignump(modulus)) { + mp_int *m = mp(modulus); + int digits = USED(m); + int bits = digits * MP_DIGIT_BIT; + int bits_needed = bits + 32; + int digits_needed = (bits_needed + MP_DIGIT_BIT - 1) / MP_DIGIT_BIT; + val out = make_bignum(); + mp_int *om = mp(out); + int i, err; + + for (i = 0; i < digits_needed; i++) { + if (i > 0) + mp_mul_2d(om, 32, om); + mp_add_d(om, rand32(r), om); + } + err = mp_mod(om, m, om); + if (err != MP_OKAY) + goto invalid; + return out; + } else if (fixnump(modulus)) { + cnum m = c_num(modulus); + if (m <= 0) + goto invalid; + return num(rand32(r) % m); + } +invalid: + uw_throwf(numeric_error_s, lit("random: invalid modulus ~s"), + modulus, nao); +} + +void rand_init(void) +{ + prot1(&random_state); + random_state_s = intern(lit("random-state"), user_package); + random_state = make_random_state(num(42)); +} |