empirical test of average rotation angle of uniform rand rotations

averot.c

// Public Domain under http://unlicense.org, see link for details.

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#include "xmmintrin.h"

//#define RANDOMIZE

// histogram size
#define HLEN 80

#define TLEN 0x7FFFFFF
// number of samples per axis of rotation
#define DLEN 512

static inline float sgn(float x) { return copysignf(1.f,x); }

// external code: xoroshiro128+

uint64_t rng_state[2];

#define TO_FP32 (1.f/16777216.f)
#define ULP1 TO_FP32
#define F32_MIN_NORMAL 1.17549435082228750796873653722224567781866555677209e-38f
#define EPS (F32_MIN_NORMAL)

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;
}
// end: xoroshiro128+

inline float rsqrtf(float v) { return 1.f/sqrtf(v); }

void reset_generators(uint64_t s0, uint64_t s1)
{
  rng_state[0] = s0;
  rng_state[1] = s1;
  rng_u64();
}

typedef struct { float x,y;     } vec2_t;
typedef union  { struct{ float x,y,z,w; }; float f[4]; } quat_t;

static inline void quat_set(quat_t* q, float x, float y, float z, float w)
{
  q->x=x; q->y=y; q->z=z; q->w=w;
}

// uniform in postive-x half disc
static float uniform_hdisc(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;
    y = 2.f*y-1.f; d += y*y;
  } while(d >= 1.f);

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

  return d;
}

static float uniform_disc(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;
}

static double uniform_quat(quat_t* q)
{
  vec2_t p0,p1;
  float  d1 = uniform_disc(&p1) + EPS;
  float  s1 = rsqrtf(d1);
  float  d0 = uniform_hdisc(&p0);
  float  s0 = sqrtf(1.f-d0);
  float  s  = s0*s1;

  quat_set(q, p0.y, s*p1.x, s*p1.y, p0.x);

  // hack to return angle of quaternion (in quaternion space)
  return acosf(p0.x);
}

int main(void)
{
  uint64_t s0;
  uint64_t s1;

#ifdef RANDOMIZE
  s0 = _rdtsc();
#else  
  s0 = 0x77535345;
#endif  
  s1 = 0x1234567;

  reset_generators(s0,s1);

  double sum = 0.0;
  quat_t q;
  uint32_t hi=0;
  uint32_t lo=0;

  // sum up the angle of TLEN random rotations.
  for(uint32_t i=0; i<TLEN; i++) {
    double a = uniform_quat(&q);

    // is angle above or below expected 'median' rotation
    if (a > 1.15494073) hi++; else lo++;

    // sum up for computing average
    sum += a;
  }

  // measure in quaternion space
  double radians = sum*(1.0/TLEN);

  // multiply by 2 to get standard angle measure
  printf("ave    = %f (radians)\n", 2.0*radians);
  printf("       = %f (radians)\n", 2.0*57.29577951308232087680*radians);
  printf("median = %f/%f\n", lo*(1.0/TLEN),hi*(1.0/TLEN));

  return 0;
}

Okay, thanks for explaining. Take it easy

No prob