uniform quaternions

gistfile1.txt

// xoroshiro128+
uint64_t rng_state[2];

#if defined(_MSC_VER)
#define TO_FP32 (1.f/16777216.f)
#else
#define TO_FP32   0x1p-24f
#endif

static inline uint64_t rotl(const uint64_t v, int i)
{
  return (v << i)|(v >> (64-i));
}

static inline uint64_t rng_u64(void)
{
  uint64_t s0 = rng_state[0];
  uint64_t s1 = rng_state[1];
  uint64_t r  = s0 + s1;
  
  s1 ^= s0;
  rng_state[0] = rotl(s0,55) ^ s1 ^ (s1<<14);
  rng_state[1] = rotl(s1,36);
  
  return r;
}

static inline float rng_f32(void)
{
  return (rng_u64() >> 40)*TO_FP32;
}

float uniform_disk(vec2_t* p)
{
  float d,x,y;
  uint64_t v;

  do {
    v = rng_u64();
    x = (v >> 40)*TO_FP32;
    y = (v & 0xFFFFFF)*TO_FP32;
    x = 2.f*x-1.f; d  = x*x;
    y = 2.f*y-1.f; d += y*y;
  } while(d >= 1.f);

  p->x = x;
  p->y = y;

  return d;
}

// uniform on quarter disk {x,y}>=0
float uniform_qdisk(vec2_t* p)
{
  float d,x,y;
  uint64_t v;

  do {
    v  = rng_u64();
    x  = (v >> 40)*TO_FP32;
    y  = (v & 0xFFFFFF)*TO_FP32;
    d  = x*x;
    d += y*y;
  } while(d >= 1.f);

  p->x = x;
  p->y = y;

  return d;
}

// uniform dist on s3
void uniform_quat(quat_t* q)
{
  // d1 really must exclude approaching (0,0)...statisticaly unlikely
  // so don't bother since this is for testing purposes. 
  vec2_t p0,p1;
  float  d0 = uniform_disk(&p0);
  float  d1 = uniform_disk(&p1);  
  float  s  = sqrtf((1.f-d0)/d1);
  quat_set(q, p0.x, p0.y, s*p1.x, s*p1.y);
}

// uniform dist with all positive components
void uniform_quat_p(quat_t* q)
{
  vec2_t p0,p1;
  float  d0 = uniform_qdisk(&p0);
  float  d1 = uniform_qdisk(&p1);
  float  s  = sqrtf((1.f-d0)/d1); // ditto 
  
  quat_set(q, p0.x, p0.y, s*p1.x, s*p1.y);
}