class_group 0.6.0

/* Copyright (C) 2000  The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
/********************************************************************/
/*                                                                  */
/*                      PSEUDO-RANDOM INTEGERS                      */
/*                                                                  */
/********************************************************************/
#include "pari.h"
#include "paripriv.h"
/********************************************************************/
/*                    XORGEN (Richard P. Brent)                     */
/*          http://wwwmaths.anu.edu.au/~brent/random.html           */
/*        (initial adaptation to PARI/GP by Randall Rathbun)        */
/********************************************************************/
/* Adapted from xorgens.c version 3.04, Richard P. Brent, 20060628 (GPL).
 * 32-bit or 64-bit integer random number generator with period at
 * least 2**4096-1. It is assumed that "ulong" is a 32-bit or 64-bit integer */

#ifdef LONG_IS_64BIT
  typedef ulong u64;
#else
  typedef unsigned long long u64;
static u64
_32to64(ulong a, ulong b) { u64 v = a; return (v<<32)|b; }
static void
_64to32(u64 v, ulong *a, ulong *b) { *a = v>>32; *b = v&0xFFFFFFFF; }
#endif
static THREAD u64 state[64];
static THREAD u64 xorgen_w;
static THREAD int xorgen_i;
/* weyl = odd approximation to 2^64*(sqrt(5)-1)/2. */
static const u64 weyl = (((u64)0x61c88646U)<<32)|((u64)0x80b583ebU);

static u64
block(void)
{
  const int r = 64;
  const int a = 33, b = 26, c = 27, d = 29, s = 53;
  u64 t, v, w;
  xorgen_i = (xorgen_i+1)&(r-1);
  t = state[xorgen_i];
  v = state[(xorgen_i+(r-s))&(r-1)];   /* index is (i-s) mod r */
  t ^= t<<a; t ^= t>>b;                   /* (I + L^a)(I + R^b) */
  v ^= v<<c; v ^= v>>d;                   /* (I + L^c)(I + R^d) */
  w = t^v;
  return state[xorgen_i] = w;               /* update circular array */
}

/* v > 0 */
static void
init_xor4096i(u64 v)
{
  const int r = 64;
  int k;

  for (k = r; k > 0; k--) {/* avoid correlations for close seeds */
    v ^= v<<10; v ^= v>>15; /* recurrence has period 2**64-1 */
    v ^= v<<4;  v ^= v>>13;
  }
  for (xorgen_w = v, k = 0; k < r; k++) { /* initialise circular array */
    v ^= v<<10; v ^= v>>15;
    v ^= v<<4;  v ^= v>>13;
    state[k] = v + (xorgen_w+=weyl);
  }
  /* discard first 4*r results */
  for (xorgen_i = r-1, k = 4*r; k > 0; k--) (void)block();
}

void pari_init_rand(void) { init_xor4096i(1UL); }

static u64
rand64(void)
{
  u64 v = block();
  xorgen_w += weyl; /* update Weyl generator */
  return v + (xorgen_w ^ (xorgen_w>>27));
}

/* One random number uniformly distributed in [0..2**BIL) is returned, where
 * BIL = 8*sizeof(ulong) = 32 or 64. */
ulong
pari_rand(void) { return rand64(); }

void
setrand(GEN x)
{
  const int r2 = numberof(state);
  long i, lx;
  u64 v;
  GEN xp;
  if (typ(x)!=t_INT) pari_err_TYPE("setrand",x);
  if (signe(x) <= 0) pari_err_DOMAIN("setrand","n", "<=", gen_0, x);
  lx = lgefint(x);
  if (lx == 3) { v = x[2]; init_xor4096i(v); return; }
#ifndef LONG_IS_64BIT
  if (lx == 4)
  {
    v = _32to64(*int_W(x,1),*int_W(x,0));
    init_xor4096i(v); return;
  }
#endif
  xp = int_LSW(x);
#ifdef LONG_IS_64BIT
  if (lx != 2 + r2+2)
    pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
  for (i = 0; i < r2; i++, xp = int_nextW(xp)) state[i] = *xp;
  xorgen_w = *xp; xp = int_nextW(xp);
#else
  if (lx != 2 + 2*r2+3)
    pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
  for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
    state[i] = _32to64(*int_nextW(xp), *xp);
  xorgen_w = _32to64(*int_nextW(xp), *xp); xp = int_nextW(int_nextW(xp));
#endif
  xorgen_i =  (*xp) & 63;
}

GEN
getrand(void)
{
  const int r2 = numberof(state);
  GEN x;
  ulong *xp;
  long i;
  if (xorgen_i < 0) init_xor4096i(1UL);

#ifdef LONG_IS_64BIT
  x = cgetipos(2+r2+2); xp = (ulong *) int_LSW(x);
  for (i = 0; i < r2; i++, xp = int_nextW(xp)) *xp = state[i];
  *xp = xorgen_w; xp = int_nextW(xp);
#else
  x = cgetipos(2+2*r2+3); xp = (ulong *) int_LSW(x);
  for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
    _64to32(state[i], int_nextW(xp), xp);
  _64to32(xorgen_w, int_nextW(xp), xp); xp = int_nextW(int_nextW(xp));
#endif
  *xp = xorgen_i? xorgen_i: 64; return x;
}

/* assume 0 <= k <= BITS_IN_LONG. Return uniform random 0 <= x < (1<<k) */
long
random_bits(long k) { return rand64() >> (64-k); }

/********************************************************************/
/*                                                                  */
/*                         GENERIC ROUTINES                         */
/*                                                                  */
/********************************************************************/

/* assume n > 0 */
ulong
random_Fl(ulong n)
{
  ulong d;
  int shift;
#ifdef LONG_IS_64BIT
  int SHIFT = 0;
#else
  int SHIFT = 32;
#endif

  if (n == 1) return 0;

  shift = bfffo(n); /* 2^(BIL-shift) > n >= 2^(BIL-shift-1)*/
  /* if N a power of 2, increment shift. No reject */
  if ((n << shift) == HIGHBIT) return rand64() >> (SHIFT+shift+1);
  for (;;) {
    d = rand64() >> (SHIFT+shift); /* d < 2^(64-shift) uniformly distributed */
    /* reject strategy: proba success = n 2^(shift-64), in [1/2, 1[ */
    if (d < n) return d;
  }
}

/* assume N > 0, see random_Fl() for algorithm. Make sure that 32-bit and
 * 64-bit architectures produce the same integers (consuming random bits
 * by packets of 64) */
GEN
randomi(GEN N)
{
  long lx = lgefint(N);
  GEN x, d;
  int shift;

  if (lx == 3) return utoi( random_Fl(N[2]) );

  shift = bfffo(*int_MSW(N));
  /* if N a power of 2, increment shift */
  if (Z_ispow2(N) && ++shift == BITS_IN_LONG) { shift = 0; lx--; }
  x = cgetipos(lx);
  for (;;) {
    GEN y, MSW = int_MSW(x), STOP = MSW;
#ifdef LONG_IS_64BIT
    for (d = int_LSW(x); d != STOP; d = int_nextW(d)) *d = rand64();
    *d = rand64() >> shift;
#else
    if (!odd(lx)) STOP = int_precW(STOP);
    /* STOP points to where MSW would in 64-bit */
    for (d = int_LSW(x); d != STOP; d = int_nextW(d))
    {
      ulong a, b; _64to32(rand64(), &a,&b);
      *d = b; d = int_nextW(d);
      *d = a;
    }
    {
      ulong a, b; _64to32(rand64() >> shift, &a,&b);
      if (d == MSW) /* 32 bits needed */
        *d = a;
      else
      { /* 64 bits needed */
        *d = b; d = int_nextW(d);
        *d = a;
      }
    }
#endif
    y = int_normalize(x, 0);
    if (abscmpii(y, N) < 0) return y;
  }
}

GEN
random_F2x(long d, long vs)
{
  ulong db, dl = dvmduBIL(d,&db);
  long i, l = 2 + dl + !!db;
  GEN y = cgetg(l,t_VECSMALL); y[1] = vs;
#ifdef LONG_IS_64BIT
  for (i=2; i<l; i++) uel(y,i) = rand64();
#else
  for (i=2; i<l-1; i+=2)
  {
    u64 v = rand64();
    uel(y,i)   = (ulong) v;
    uel(y,i+1) = (ulong) (v>>32);
  }
  if (i<l) uel(y,i) = (ulong) rand64();
#endif
  if (db) uel(y,l-1) &= ((1UL<<db)-1UL);
  return F2x_renormalize(y,l);
}

GEN
randomr(long prec)
{
  pari_sp av;
  long b;
  GEN x, y;
  if (prec <= 2) return real_0_bit(0);
  x = cgetr(prec); av = avma;
  b = prec2nbits(prec);
  y = randomi(int2n(b));
  if (!signe(y)) return real_0_bit(b);
  affir(y, x); shiftr_inplace(x, - b);
  avma = av; return x;
}

static GEN
polrandom(GEN N) /* assume N!=0 */
{
  long i, d = lg(N);
  GEN z = leading_coeff(N);
  GEN y = cgetg(d,t_POL);
  y[1] = evalsigne(1) | evalvarn(varn(N));
  for (i=2; i<d; i++) gel(y,i) = genrand(z);
  return normalizepol_lg(y,d);
}

GEN
genrand(GEN N)
{
  GEN z;
  if (!N) return utoi( random_bits(31) );
  switch(typ(N))
  {
    case t_INT:
      if (signe(N)<=0) pari_err_DOMAIN("random","N","<=",gen_0,gen_0);
      return randomi(N);
    case t_REAL:
      return randomr(realprec(N));
    case t_INTMOD:
      z = cgetg(3, t_INTMOD);
      gel(z,1) = icopy(gel(N,1));
      gel(z,2) = randomi(gel(N,1)); return z;
    case t_FFELT:
      return ffrandom(N);
    case t_POL:
      if (signe(N)==0) return pol_0(varn(N));
      return polrandom(N);
    case t_VEC:
      if (lg(N) == 3)
      {
        pari_sp av = avma;
        GEN a = gel(N,1), b = gel(N,2), d;
        if (typ(a) != t_INT) a = gceil(a);
        if (typ(b) != t_INT) b = gfloor(b);
        if (typ(a) != t_INT || typ(b) != t_INT) pari_err_TYPE("random", N);
        d = subii(b,a);
        if (signe(d) < 0) pari_err_TYPE("random([a,b]) (a > b)", N);
        return gerepileuptoint(av, addii(a, randomi(addiu(d,1))));
      }
      return ellrandom(N);
    default:
      pari_err_TYPE("genrand",N);
      return NULL;/*LCOV_EXCL_LINE*/
  }
}