Efficient generation of chi random variables I need to generate random variables generated from a chi distribution (not chi-squared!).  There doesn't seem to be standard mechanism in C++ in (for example) Boost::Random and hence I am looking for an alternative implementation implemented in C++ (or an easily compatible language such as C or FORTRAN).
There is one paper in ACM TOMS on the topic. Or, should I simply take the square root of a chi-squared random variable, of which there are several implementations?
Which is likely to be fastest - unfortunately, I need to draw a lot of them.

Update:
In response to the below comments, these are some notes for my specific application:


*

*I am looking for a C/C++ compatible implementation, or sufficient detail to program one myself.

*The application is CPU sensitive, and is run on a modest embedded platform (400MHz with a scalar double FPU and limited cache, but no sqrt(), sin(), exp() or similar instructions). Since actual run-time is not easy to predict a priori, I am looking for algorithms that are amongst the fastest and known to perform well. I can evaluate each of these on my own hardware.

*I already have facilities for generating uniform and normal random numbers. I can get chi-squared if necessary.

*DoF is always an integer, usually between 15 and 39.


The answer I accept will have the following content:


*

*Which methods are considered "fast", and how are they implemented?

*Do they otherwise have known disadvantages (e.g. numerical stability)?

*Are there any methods that should be specifically avoided?

 A: Here's a pure C implementation that generates a little over three million samples per second (using df=20) on my machine (Macbook Pro i7 @ 2.7 GHz).
I've assumed the only mathematical functions supported natively are multiplication, division, addition, subtraction, and the rand() function from the standard library. You can easily replace the rand() function with something equivalent - it seems to be the bottleneck in this program after profiling. I've avoided linking to math.h to avoid pulling in the extra code.
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
const int HALF_RAND_MAX = RAND_MAX/2;
const float c1 = -3.0/2;
const float c2 = 11.0/6;
const float c3 = -25.0/12;
const float c4 = 137.0/60;
const float c5 = -49.0/20;

// Reuse of the other squared random normal.
static int hasSpare = 0;
static float spare = 0.0;

float randu() {
  // Generates a random float uniformly between [-1.0, 1.0) using RAND_MAX
  return ((float)rand())/(HALF_RAND_MAX) - 1;
}

inline float lnxinvx(float x) {
  // return approximation for ln(x) / x
  float v = x - 1;
  return v + v*(c1 + v*(c2 + v*(c3 + v*(c4 + v*c5))));
}
float fast_sqrt(float x) {
  // Adapted from Carmac's famous Q_sqrt trick.
  int i;
  float xhalf = x*0.5f;
  i = *(int*)&x;
  i = 0x5f3759df - (i >> 1);
  x = *(float*)&i;
  x = x*(1.5f - (xhalf * x * x));
  x = x*(1.5f - (xhalf * x * x));
  return 1.0f/x;
}

float squared_random_normal() {
  // Generate a squared random normal variable, or use the spare one.
  if (hasSpare) {
    hasSpare = 0;
    return spare;
  } else {
    float u = randu();
    float v = randu();
    float s = u*u + v*v;
    while (s >= 1) {
      u = randu();
      v = randu();
      s = u*u + v*v;
    }
    float mul = -2.0*lnxinvx(s);
    spare = v * v * mul;
    hasSpare = 1;
    return u * u * mul;
  }
}

float rand_chi(unsigned int df) {
  // Add up "df" squared random normals and take the square root.
  float acc = 0;
  int i = 0;
  for(i; i < df; ++i) {
    acc += squared_random_normal();
  }
  return fast_sqrt(acc);
}

int main(char* argc, char** argv) {
  srand(time(NULL));
  int i;
  int df = atoi(argv[1]);
  int n = atoi(argv[2]);
  printf("Using %i degrees of freedom\n", df);
  printf("Generating %i samples\n", n);
  for(i = 0; i < n; ++i) {
    //printf("Next variate: %f\n", rand_chi(df));
    rand_chi(df);
  }
}

And here's the timing information, for 20 degrees of freedom.
Wed Aug 14 00:08:04 ~/Dropbox/RandomCodeSnippets/chidist $ gcc -O3 -funroll-loops -fwhole-program chidist.c -o chidist
Wed Aug 14 00:08:09 ~/Dropbox/RandomCodeSnippets/chidist $ time ./chidist 20 10000000
Using 20 degrees of freedom
Generating 10000000 samples

real    0m3.208s
user    0m3.205s
sys     0m0.002s

Hopefully this is roughly what you were looking for (or is roughly fast enough).
Thanks!
A: To generate Chi distribution random variate in Mathematica, general expression is
RandomVariate[ChiDistribution[degree of freedom], number of random variates]

For example we have to generate five random variates with three degree of freedom the code is
RandomVariate[ChiDistribution[3],5]

The output is
{0.753102, 1.6647, 1.25129, 0.456877, 0.632508}

A: I haven't yet used RcppArmadillo package but you may look into it. This package is useful if C++ has been decided as the language of choice. 
