chiark - git - mdw - catacomb/blob - math/mprand.c

   1 /* -*-c-*-
   2  *
   3  * Generate a random multiprecision integer
   4  *
   5  * (c) 1999 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include <mLib/alloc.h>
  31
  32 #include "grand.h"
  33 #include "mp.h"
  34 #include "mprand.h"
  35
  36 /*----- Main code ---------------------------------------------------------*/
  37
  38 /* --- @mprand@ --- *
  39  *
  40  * Arguments:   @mp *d@ = destination integer
  41  *              @unsigned b@ = number of bits
  42  *              @grand *r@ = pointer to random number source
  43  *              @mpw or@ = mask to OR with low-order bits
  44  *
  45  * Returns:     A random integer with the requested number of bits.
  46  *
  47  * Use:         Constructs an arbitrarily large pseudorandom integer.
  48  *              Assuming that the generator @r@ is good, the result is
  49  *              uniformly distributed in the interval %$[2^{b - 1}, 2^b)$%.
  50  *              The result is then ORred with the given @or@ value.  This
  51  *              will often be 1, to make the result odd.
  52  *
  53  *              The length @b@ may be zero; but %$\texttt{or} \ge 2^b$% is
  54  *              not permitted.
  55  */
  56
  57 mp *mprand(mp *d, unsigned b, grand *r, mpw or)
  58 {
  59   size_t sz = (b + 7) >> 3;
  60   arena *a = (d && (d->f & MP_BURN)) ? arena_secure : arena_global;
  61   octet *v;
  62   unsigned m;
  63
  64   assert(b >= MPW_BITS || !(or >> b));
  65
  66   /* --- Special case --- */
  67
  68   if (!b) return (MP_ZERO);
  69
  70   /* --- Fill buffer with random data --- */
  71
  72   v = x_alloc(a, sz);
  73   r->ops->fill(r, v, sz);
  74
  75   /* --- Force into the correct range --- *
  76    *
  77    * This is slightly tricky.  Oh, well.
  78    */
  79
  80   b = (b - 1) & 7;
  81   m = (1 << b);
  82   v[0] = (v[0] & (m - 1)) | m;
  83
  84   /* --- Mask, load and return --- */
  85
  86   d = mp_loadb(d, v, sz);
  87   if (or) {
  88     assert(d->sz);
  89     if (!MP_LEN(d)) d->vl = d->v + 1;
  90     d->v[0] |= or;
  91   }
  92   memset(v, 0, sz);
  93   x_free(a, v);
  94   return (d);
  95 }
  96
  97 /* --- @mprand_range@ --- *
  98  *
  99  * Arguments:   @mp *d@ = destination integer
 100  *              @mp *l@ = limit for random number
 101  *              @grand *r@ = random number source
 102  *              @mpw or@ = mask for low-order bits
 103  *
 104  * Returns:     A pseudorandom integer, unformly distributed over the
 105  *              interval %$[0, l)$%.
 106  *
 107  * Use:         Generates a uniformly-distributed pseudorandom number in the
 108  *              appropriate range.  We must have %$l > 0$%.
 109  */
 110
 111 mp *mprand_range(mp *d, mp *l, grand *r, mpw or)
 112 {
 113   size_t b = mp_bits(l);
 114   size_t sz = (b + 7) >> 3;
 115   arena *a = (d && (d->f & MP_BURN)) ? arena_secure : arena_global;
 116   octet *v = x_alloc(a, sz);
 117   unsigned m;
 118
 119   /* --- The algorithm --- *
 120    *
 121    * Rather simpler than most.  Find the number of bits in the number %$l$%
 122    * (i.e., the integer %$b$% such that %$2^{b - 1} \le l < 2^b$%), and
 123    * generate pseudorandom integers with %$n$% bits (but not, unlike in the
 124    * function above, with the top bit forced to 1).  If the integer is
 125    * greater than or equal to %$l$%, try again.
 126    *
 127    * This is similar to the algorithms used in @lcrand_range@ and friends,
 128    * except that I've forced the `raw' range of the random numbers such that
 129    * %$l$% itself is the largest multiple of %$l$% in the range (since, by
 130    * the inequality above, %$2^b \le 2l$%).  This removes the need for costly
 131    * division and remainder operations.
 132    *
 133    * As usual, the number of iterations expected is two.
 134    */
 135
 136   assert(MP_POSP(l));
 137   b = ((b - 1) & 7) + 1;
 138   m = (1 << b) - 1;
 139   do {
 140     r->ops->fill(r, v, sz);
 141     v[0] &= m;
 142     d = mp_loadb(d, v, sz);
 143     d->v[0] |= or;
 144   } while (MP_CMP(d, >=, l));
 145
 146   /* --- Done --- */
 147
 148   memset(v, 0, sz);
 149   x_free(a, v);
 150   return (d);
 151 }
 152
 153 /*----- That's all, folks -------------------------------------------------*/