rand_custom.c

/*
rand_custom.c is a rewritten version of the glibc/stdlib/random_r.c in order to get the system state of seeding threadprivate
Copyright (C) 2017  Mollenhauer, Ruhl

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library 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
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, see
<http://www.gnu.org/licenses/>.
*/

/*
 * This is derived from this source:
 * https://github.com/lattera/glibc/blob/a2f34833b1042d5d8eeb263b4cf4caaea138c4ad/stdlib/random_r.c
 * it is rewritten, in order to get the system state of the random seed threadprivate, by Mollenhauer & Ruhl
 */

#include "rand_custom.h"

static struct random_data unsafe_state = {0};

int __srandom_r_custom(unsigned int seed, struct random_data *buf) {
    int type;
    int32_t *state;
    long int i;
    int32_t word;
    int32_t *dst;
    int kc;

    if (buf == NULL)
        goto fail;
    type = buf->rand_type;
    if ((unsigned int) type >= MAX_TYPES)
        goto fail;

    state = buf->state;
    /* We must make sure the seed is not 0.  Take arbitrarily 1 in this case.  */
    if (seed == 0)
        seed = 1;
    state[0] = seed;
    if (type == TYPE_0)
        goto done;

    dst = state;
    word = seed;
    kc = buf->rand_deg;
    for (i = 1; i < kc; ++i) {
        /* This does:
         state[i] = (16807 * state[i - 1]) % 2147483647;
       but avoids overflowing 31 bits.  */
        long int hi = word / 127773;
        long int lo = word % 127773;
        word = 16807 * lo - 2836 * hi;
        if (word < 0)
            word += 2147483647;
        *++dst = word;
    }

    buf->fptr = &state[buf->rand_sep];
    buf->rptr = &state[0];
    kc *= 10;
    while (--kc >= 0) {
        int32_t discard;
        (void) __random_r_custom(buf, &discard);
    }

    done:
    return 0;

    fail:
    return -1;
}

void srand_custom(unsigned int x) {
    static int32_t randtbl[DEG_3 + 1] =
            {
                    TYPE_3,

                    -1726662223, 379960547, 1735697613, 1040273694, 1313901226,
                    1627687941, -179304937, -2073333483, 1780058412, -1989503057,
                    -615974602, 344556628, 939512070, -1249116260, 1507946756,
                    -812545463, 154635395, 1388815473, -1926676823, 525320961,
                    -1009028674, 968117788, -123449607, 1284210865, 435012392,
                    -2017506339, -911064859, -370259173, 1132637927, 1398500161,
                    -205601318,
            };

    // Set internal state and random table to be thread private
#pragma omp threadprivate(randtbl)
#pragma omp threadprivate(unsafe_state)

    unsafe_state.fptr = &randtbl[SEP_3 + 1];
    unsafe_state.fptr = &randtbl[1];
    unsafe_state.state = &randtbl[1];
    unsafe_state.rand_type = TYPE_3;
    unsafe_state.rand_deg = DEG_3;
    unsafe_state.rand_sep = SEP_3;
    unsafe_state.end_ptr = &randtbl[sizeof(randtbl) / sizeof(randtbl[0])];
    (void) __srandom_r_custom(x, &unsafe_state);
}

int __random_r_custom(struct random_data *buf, int32_t *result) {
    int32_t *state;

    if (buf == NULL || result == NULL)
        goto fail;

    state = buf->state;

    if (buf->rand_type == TYPE_0) {
        int32_t val = state[0];
        val = ((state[0] * 1103515245) + 12345) & 0x7fffffff;
        state[0] = val;
        *result = val;
    } else {
        int32_t *fptr = buf->fptr;
        int32_t *rptr = buf->rptr;
        int32_t *end_ptr = buf->end_ptr;
        int32_t val;

        val = *fptr += *rptr;
        /* Chucking least random bit.  */
        *result = (val >> 1) & 0x7fffffff;
        ++fptr;
        if (fptr >= end_ptr) {
            fptr = state;
            ++rptr;
        } else {
            ++rptr;
            if (rptr >= end_ptr)
                rptr = state;
        }
        buf->fptr = fptr;
        buf->rptr = rptr;
    }
    return 0;

    fail:
    errno = EINVAL;
    return -1;
}

int32_t rand_custom() {
    int32_t retval;
    (void) __random_r_custom(&unsafe_state, &retval);
    return retval;
}

rand_custom.h

#pragma once

#include <stdio.h>
#include <stdint.h>
#include <asm/errno.h>
#include <errno.h>

/* Linear congruential.  */
#define    TYPE_0        0
#define    BREAK_0        8
#define    DEG_0        0
#define    SEP_0        0

/* x**7 + x**3 + 1.  */
#define    TYPE_1        1
#define    BREAK_1        32
#define    DEG_1        7
#define    SEP_1        3

/* x**15 + x + 1.  */
#define    TYPE_2        2
#define    BREAK_2        64
#define    DEG_2        15
#define    SEP_2        1

/* x**31 + x**3 + 1.  */
#define    TYPE_3        3
#define    BREAK_3        128
#define    DEG_3        31
#define    SEP_3        3

/* x**63 + x + 1.  */
#define    TYPE_4        4
#define    BREAK_4        256
#define    DEG_4        63
#define    SEP_4    1


struct random_data
{
    int32_t *fptr;		/* Front pointer.  */
    int32_t *rptr;		/* Rear pointer.  */
    int32_t *end_ptr;		/* Pointer behind state table.  */
    int32_t *state;
    int rand_deg;
    int rand_sep;
    int rand_type;
};

/* Array versions of the above information to make code run faster.
   Relies on fact that TYPE_i == i.  */

#define    MAX_TYPES    5    /* Max number of types above.  */

int __random_r_custom(struct random_data *buf, int32_t *result);

int __srandom_r_custom(unsigned int seed, struct random_data *buf);

void srand_custom(unsigned int x);

int32_t rand_custom();