Merge branch '2.4.x' into 2.5.x

[catacomb] / math / mprand.c

/* -*-c-*-
 *
 * Generate a random multiprecision integer
 *
 * (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 <mLib/alloc.h>

#include "grand.h"
#include "mp.h"
#include "mprand.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @mprand@ --- *
 *
 * Arguments:	@mp *d@ = destination integer
 *		@unsigned b@ = number of bits
 *		@grand *r@ = pointer to random number source
 *		@mpw or@ = mask to OR with low-order bits
 *
 * Returns:	A random integer with the requested number of bits.
 *
 * Use:		Constructs an arbitrarily large pseudorandom integer.
 *		Assuming that the generator @r@ is good, the result is
 *		uniformly distributed in the interval %$[2^{b - 1}, 2^b)$%.
 *		The result is then ORred with the given @or@ value.  This
 *		will often be 1, to make the result odd.
 *
 *		The length @b@ may be zero; but %$\texttt{or} \ge 2^b$% is
 *		not permitted.
 */

mp *mprand(mp *d, unsigned b, grand *r, mpw or)
{
  size_t sz = (b + 7) >> 3;
  arena *a = (d && (d->f & MP_BURN)) ? arena_secure : arena_global;
  octet *v;
  unsigned m;

  assert(b >= MPW_BITS || !(or >> b));

  /* --- Special case --- */

  if (!b) return (MP_ZERO);

  /* --- Fill buffer with random data --- */

  v = x_alloc(a, sz);
  r->ops->fill(r, v, sz);

  /* --- Force into the correct range --- *
   *
   * This is slightly tricky.  Oh, well.
   */

  b = (b - 1) & 7;
  m = (1 << b);
  v[0] = (v[0] & (m - 1)) | m;

  /* --- Mask, load and return --- */

  d = mp_loadb(d, v, sz);
  if (or) {
    assert(d->sz);
    if (!MP_LEN(d)) d->vl = d->v + 1;
    d->v[0] |= or;
  }
  memset(v, 0, sz);
  x_free(a, v);
  return (d);
}

/* --- @mprand_range@ --- *
 *
 * Arguments:	@mp *d@ = destination integer
 *		@mp *l@ = limit for random number
 *		@grand *r@ = random number source
 *		@mpw or@ = mask for low-order bits
 *
 * Returns:	A pseudorandom integer, unformly distributed over the
 *		interval %$[0, l)$%.
 *
 * Use:		Generates a uniformly-distributed pseudorandom number in the
 *		appropriate range.  We must have %$l > 0$%.
 */

mp *mprand_range(mp *d, mp *l, grand *r, mpw or)
{
  size_t b = mp_bits(l);
  size_t sz = (b + 7) >> 3;
  arena *a = (d && (d->f & MP_BURN)) ? arena_secure : arena_global;
  octet *v = x_alloc(a, sz);
  unsigned m;

  /* --- The algorithm --- *
   *
   * Rather simpler than most.  Find the number of bits in the number %$l$%
   * (i.e., the integer %$b$% such that %$2^{b - 1} \le l < 2^b$%), and
   * generate pseudorandom integers with %$n$% bits (but not, unlike in the
   * function above, with the top bit forced to 1).  If the integer is
   * greater than or equal to %$l$%, try again.
   *
   * This is similar to the algorithms used in @lcrand_range@ and friends,
   * except that I've forced the `raw' range of the random numbers such that
   * %$l$% itself is the largest multiple of %$l$% in the range (since, by
   * the inequality above, %$2^b \le 2l$%).  This removes the need for costly
   * division and remainder operations.
   *
   * As usual, the number of iterations expected is two.
   */

  assert(MP_POSP(l));
  b = ((b - 1) & 7) + 1;
  m = (1 << b) - 1;
  do {
    r->ops->fill(r, v, sz);
    v[0] &= m;
    d = mp_loadb(d, v, sz);
    d->v[0] |= or;
  } while (MP_CMP(d, >=, l));

  /* --- Done --- */

  memset(v, 0, sz);
  x_free(a, v);
  return (d);
}

/*----- That's all, folks -------------------------------------------------*/