contrib/mpfr/random2.c

   1 /* mpfr_random2 -- Generate a positive random mpfr_t of specified size, with
   2    long runs of consecutive ones and zeros in the binary representation.
   3    Intended for testing.
   4
   5 Copyright 1999, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   6 Contributed by the Arenaire and Cacao projects, INRIA.
   7
   8 This file is part of the GNU MPFR Library.
   9
  10 The GNU MPFR Library is free software; you can redistribute it and/or modify
  11 it under the terms of the GNU Lesser General Public License as published by
  12 the Free Software Foundation; either version 2.1 of the License, or (at your
  13 option) any later version.
  14
  15 The GNU MPFR Library is distributed in the hope that it will be useful, but
  16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  18 License for more details.
  19
  20 You should have received a copy of the GNU Lesser General Public License
  21 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  22 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  23 MA 02110-1301, USA. */
  24
  25 #define MPFR_NEED_LONGLONG_H
  26 #include "mpfr-impl.h"
  27
  28 #define LOGBITS_PER_BLOCK 4
  29 #if GMP_NUMB_BITS < 32
  30 #define BITS_PER_RANDCALL GMP_NUMB_BITS
  31 #else
  32 #define BITS_PER_RANDCALL 32
  33 #endif
  34
  35 static void
  36 mpfr_random2_raw (mpfr_ptr x, mp_size_t size, mp_exp_t exp,
  37                   gmp_randstate_t rstate)
  38 {
  39   mp_size_t xn, k, ri;
  40   unsigned long sh;
  41   mp_ptr xp;
  42   mp_limb_t elimb, ran, acc;
  43   int ran_nbits, bit_pos, nb;
  44
  45   MPFR_CLEAR_FLAGS (x);
  46
  47   if (MPFR_UNLIKELY(size == 0))
  48     {
  49       MPFR_SET_ZERO (x);
  50       MPFR_SET_POS (x);
  51       return ;
  52     }
  53   else if (size > 0)
  54     {
  55       MPFR_SET_POS (x);
  56     }
  57   else
  58     {
  59       MPFR_SET_NEG (x);
  60       size = -size;
  61     }
  62
  63   xn = MPFR_LIMB_SIZE (x);
  64   xp = MPFR_MANT (x);
  65   if (size > xn)
  66     size = xn;
  67   k = xn - size;
  68
  69   /* Code extracted from GMP, function mpn_random2, to avoid the use
  70      of GMP's internal random state in MPFR */
  71
  72   _gmp_rand (&elimb, rstate, BITS_PER_RANDCALL);
  73   ran = elimb;
  74
  75   /* Start off at a random bit position in the most significant limb.  */
  76   bit_pos = GMP_NUMB_BITS - 1;
  77   ran >>= 6;                            /* Ideally   log2(GMP_NUMB_BITS) */
  78   ran_nbits = BITS_PER_RANDCALL - 6;    /* Ideally - log2(GMP_NUMB_BITS) */
  79
  80   /* Bit 0 of ran chooses string of ones/string of zeroes.
  81      Make most significant limb be non-zero by setting bit 0 of RAN.  */
  82   ran |= 1;
  83   ri = xn - 1;
  84
  85   acc = 0;
  86   while (ri >= k)
  87     {
  88       if (ran_nbits < LOGBITS_PER_BLOCK + 1)
  89         {
  90           _gmp_rand (&elimb, rstate, BITS_PER_RANDCALL);
  91           ran = elimb;
  92           ran_nbits = BITS_PER_RANDCALL;
  93         }
  94
  95       nb = (ran >> 1) % (1 << LOGBITS_PER_BLOCK) + 1;
  96       if ((ran & 1) != 0)
  97         {
  98           /* Generate a string of nb ones.  */
  99           if (nb > bit_pos)
 100             {
 101               xp[ri--] = acc | (((mp_limb_t) 2 << bit_pos) - 1);
 102               bit_pos += GMP_NUMB_BITS;
 103               bit_pos -= nb;
 104               acc = ((~(mp_limb_t) 1) << bit_pos) & GMP_NUMB_MASK;
 105             }
 106           else
 107             {
 108               bit_pos -= nb;
 109               acc |= (((mp_limb_t) 2 << nb) - 2) << bit_pos;
 110             }
 111         }
 112       else
 113         {
 114           /* Generate a string of nb zeroes.  */
 115           if (nb > bit_pos)
 116             {
 117               xp[ri--] = acc;
 118               acc = 0;
 119               bit_pos += GMP_NUMB_BITS;
 120             }
 121           bit_pos -= nb;
 122         }
 123       ran_nbits -= LOGBITS_PER_BLOCK + 1;
 124       ran >>= LOGBITS_PER_BLOCK + 1;
 125     }
 126
 127   /* Set mandatory most significant bit.  */
 128   /* xp[xn - 1] |= MPFR_LIMB_HIGHBIT; */
 129
 130   if (k != 0)
 131     {
 132       /* Clear last limbs */
 133       MPN_ZERO (xp, k);
 134     }
 135   else
 136     {
 137       /* Mask off non significant bits in the low limb.  */
 138       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (x));
 139       xp[0] &= ~MPFR_LIMB_MASK (sh);
 140     }
 141
 142   /* Generate random exponent.  */
 143   _gmp_rand (&elimb, RANDS, BITS_PER_MP_LIMB);
 144   exp = ABS (exp);
 145   MPFR_SET_EXP (x, elimb % (2 * exp + 1) - exp);
 146
 147   return ;
 148 }
 149
 150 void
 151 mpfr_random2 (mpfr_ptr x, mp_size_t size, mp_exp_t exp)
 152 {
 153   mpfr_random2_raw (x, size, exp, RANDS);
 154 }