sketch for part of constant-time random double algorithm.

gistfile1.c

#include <stdio.h>
#include <stdint.h>
/*
 * "See https://github.com/zackw/stegotorus/blob/master/src/rng.cc#L86 and the
 * paper at http://allendowney.com/research/rand/ " -- zackw.
 *
 */

/* Here's what Zack's code does: */
#if 0
  /* This is the core of Downey's algorithm: 50% of the time we
     should generate the highest exponent of a number in (0,1) (note
     that _neither_ endpoint is included right now).  25% of the
     time, we should generate the second highest exponent, 12.5% of
     the time, we should generate the third highest, and so on.  In
     other words, we should start with the highest exponent, flip a
     coin, and keep subtracting 1 until either we hit zero or the
     coin comes up heads.

     If anyone knows how to do this in _constant_ time, instead of
     variable time bounded by a constant, please tell me.
  */

  uint32_t exponent = 0x3FE; /* 1111111110 = 2^{-1} */
  do {
    if (bits.get()) break;
  } while (--exponent);


  /* Finally a slight adjustment: if the mantissa is zero, then
     half the time we should increment the exponent by one.
     Do this unconditionally if the exponent is also zero
     (so we never generate 0.0). */
  if (mantissa == 0 && (exponent == 0 || bits.get()))
    exponent++;
#endif

/* Okay, here's the start of a not-quite-right constant-time version. */

/* Compute the log2 of n.  The highest value is 63; log2_64(0) returns 0. */
unsigned
log2_64(uint64_t v)
{
  /* XXXX need to check that this _is_ constant-time! Some compilers might
   * do some silly stuff here when they see the > operations.
   *
   * Source: http://graphics.stanford.edu/~seander/bithacks.html
   */
  uint64_t r;
  unsigned int shift;

  r =     (v > 0xFFFFFFFF) << 8; v >>= r;
  shift = (v > 0xFFFF    ) << 4; v >>= shift; r |= shift;
  shift = (v > 0xFF      ) << 3; v >>= shift; r |= shift;
  shift = (v > 0xF       ) << 2; v >>= shift; r |= shift;
  shift = (v > 0x3       ) << 1; v >>= shift; r |= shift; r |= (v >> 1);

  return (unsigned)r;
}

/* Requires: u is less than 2^24.
 *
 * Returns: 1 if u is 0, 0 otherwise. */
unsigned
small_unsigned_is_zero(unsigned u)
{
  return ((u-1)>>24) & 1;
}

/* Returns 1 if u is 0, 0 otherwise. */
unsigned
uint64_is_zero(uint64_t u)
{
  /* t is less than 2^32, and has bits set iff u has any bits set. */
  const uint64_t t = (u & 0xffffffff) | (u >> 32);

  return (unsigned)( ((t-1)>>32) & 1);
}

/* Given an array of uint64_t stored in big-endian order, return the log2 of
 * the whole array. */
int
log2_many(uint64_t *array, int array_len)
{
  unsigned r = 0;
  int i;

  for (i = array_len-1; i >= 0; --i) {
    /* we want to do:

       if (array[i])
         r = log2(array[i]) + i*64;

      So here's how to do that in constant-time.
    */

    unsigned log2 = log2_64(array[i]) + (i<<6);
    unsigned elt_is_zero = uint64_is_zero(array[i]);
    /* mask== 0 if array[i]==0,
     * mask==~0    otherwise. */
    unsigned mask = elt_is_zero - 1;

    /* Again, verify that this formulation is constant-time on your compiler */
    r = (mask & r) | (~mask & log2);
  }

  return r;
}


#if 0
{

  /* Set this parameter to N such that events occurring with P< 2^(-64 * N)
   * are "impossible".  I don't believe in P==2^-128 coincidences, so I picked
   * N==2.  You might pick 1 if you're brave, or 4 if you're a bit nutty,
   * or 16 if you simply don't believe in probability.
   *
   * (If you pick PRECISION==16, this algorithm won't quite work, since it's
   * possible in theory to get a negative exponent if you do that.  I wouldn't
   * worry about that, since PRECISION==16 is for weirdos IMO.)
   **/
#define PRECISION 2

  uint64_t bits[precision];
  randomize(bits); /* through whatever uniform random bit generator we have. */

  /* We're going to take the base-2 logarithm of "bits", considering "bits" as
   * a big-endian integer. */
  unsigned log2 = log2_many(bits, precision);

#define EXP_MAX 0x3FE

  /* Now, log2 == (PRECISION*64 -1) with probability == 1/2,
          log2 == (PRECISION*64 -2) with probability == 1/4,
          log2 == (PRECISION*64 -3) with probability == 1/8....

      We want exponent to equal EXP_MAX   with probability 1/2,
                                EXP_MAX-1 with probability 1/4,
                                EXP_MAX-2 with probability 1/8...
  */
#define LOG2_MAX (PRECISION * 64 - 1)

  unsigned exponent = log2 + EXP_MAX - PRECISION_MAX;

  /* We could use the fact that mantissa < 2^52 to optimize this a bit */
  unsigned mantissa_is_zero = uint64_is_zero(mantissa);

  /* The possibility of exponent==0 is negligible (if you run the full
     algorithm) and zero (if you use the PRECISION<16 cheat above), so
     you could simplify this by removing the check on exponent..
  */
  exponent +=
    mantissa_is_zero & (randombit() | small_unsigned_is_zero(exponent));
}
#endif