struct stable_random_state {
  unsigned shift:  32;
  unsigned carry:   8;
  unsigned current: 8;
  unsigned prev:    8;
  unsigned linear:  8;
};

#define STABLE_RANDOM_NEXT_LINEAR(s) ((s) -> linear *= 73, (s) -> linear += 29, (s) -> linear)

unsigned char stable_random (struct stable_random_state * state) {
  unsigned p;
  const unsigned char cycle_start_points[] = {1, 2, 4, 8, 13, 17, 23, 26, 29, 58, 0};
  const unsigned char short_cycles[] = {0x72, 0x4f, 0x9f, 0x7b, 0x1a, 0x7b, 0x84, 0xe5, 0x56, 0x8d, 0xb0, 0x32, 0, 0, 1};
  if (!state -> shift) for (p = 0; p < 4; p ++) state -> shift = (state -> shift << 8) | STABLE_RANDOM_NEXT_LINEAR(state);
  state -> shift ^= state -> shift >> 8;
  state -> shift ^= state -> shift << 9;
  state -> shift ^= state -> shift >> 23;
  if (state -> prev || state -> current)
    for (p = 0; p < (sizeof short_cycles - 3); p += 3) {
      if ((state -> prev == short_cycles[p]) && (state -> current == short_cycles[p + 1]) && (state -> carry == short_cycles[p + 2])) {
        state -> prev = short_cycles[p + 3];
        state -> current = short_cycles[p + 4];
        state -> carry = short_cycles[p + 5];
        break;
      }
    }
  else
    for (p = 0; p < (sizeof cycle_start_points - 1); p ++) if (state -> carry == cycle_start_points[p]) {
      state -> carry = cycle_start_points[p + 1];
      if (!state -> carry) {
        state -> prev = *short_cycles;
        state -> current = short_cycles[1];
        state -> carry = short_cycles[2];
        STABLE_RANDOM_NEXT_LINEAR(state);
      }
      break;
    }
  if (state -> carry >= 210) state -> carry -= 210;
  p = state -> carry + state -> prev + state -> current;
  if (!p || (p == 719)) {
    state -> prev = STABLE_RANDOM_NEXT_LINEAR(state);
    state -> carry = STABLE_RANDOM_NEXT_LINEAR(state);
    state -> current = STABLE_RANDOM_NEXT_LINEAR(state);
    if (state -> carry >= 210) state -> carry -= 210;
  }
  p = 210 * state -> prev + state -> carry;
  state -> prev = state -> current;
  state -> current = p & 0xff;
  state -> carry = p >> 8;
  STABLE_RANDOM_NEXT_LINEAR(state);
  p = state -> shift >> ((state -> linear >> 3) & 24);
  switch ((state -> linear >> 4) & 3) {
    case 0:
      return p + state -> current;
    case 1:
      return p ^ state -> current;
    case 2:
      return p - state -> current;
    default:
      return state -> current - p;
  }
}

void stable_random_initialize (struct stable_random_state * state, unsigned seed1, unsigned seed2, unsigned initial_mixing_count) {
  // assuming here that unsigned is 32-bit
  state -> shift = seed1;
  1[(unsigned *) state] = seed2;
  while (initial_mixing_count --) stable_random(state);
}