```/****************************************************************************/
/* MATHFNS.C: utility routines for various sorts of math operations. Most   */
/* these functions work with real values, meaning that they can handle      */
/* either floats or doubles, depending on compiler switches.                */
/* Copyright (c) 1999 by Joshua E. Barnes, Tokyo, JAPAN.                    */
/****************************************************************************/

#include "stdinc.h"
#include "mathfns.h"

long random(void);

/*
* RSQR, RQBE: compute x*x and x*x*x.
*/

real rsqr(real x)
{
return (x * x);
}

real rqbe(real x)
{
return (x * x * x);
}

/*
* RLOG2, REXP2: log, inverse log to base two.
*/

real rlog2(real x)
{
return (rlog(x) / M_LN2);
}

real rexp2(real x)
{
return (rexp(M_LN2 * x));
}

/*
* RDEX: inverse log base ten.
*/

real rdex(real x)
{
return (rexp(M_LN10 * x));
}

#if defined(SINGLEPREC)

/*
* FCBRT: floating cube root.
*/

float fcbrt(float x)
{
return ((float) cbrt((double) x));
}

#endif

/*
* XRANDOM: floating-point random number routine.
*/

double xrandom(double xl, double xh)
{

return (xl + (xh - xl) * ((double) random()) / 2147483647.0);
}

/*
* GRANDOM: normally distributed random number (polar method).
* Reference: Knuth, vol. 2, p. 104.
*/

double grandom(double mean, double sdev)
{
double v1, v2, s;

do {
v1 = xrandom(-1.0, 1.0);
v2 = xrandom(-1.0, 1.0);
s = v1*v1 + v2*v2;
} while (s >= 1.0);
return (mean + sdev * v1 * sqrt(-2.0 * log(s) / s));
}

/*
* PICKSHELL: pick point on shell.
*/

void pickshell(real vec[], int ndim, real rad)
{
real rsq, rscale;
int i;

do {
rsq = 0.0;
for (i = 0; i < ndim; i++) {
vec[i] = xrandom(-1.0, 1.0);
rsq = rsq + vec[i] * vec[i];
}
} while (rsq > 1.0);
for (i = 0; i < ndim; i++)
vec[i] = vec[i] * rscale;
}

/*
* PICKBALL: pick point within ball.
*/

void pickball(real vec[], int ndim, real rad)
{
real rsq;
int i;

do {
rsq = 0.0;
for (i = 0; i < ndim; i++) {
vec[i] = xrandom(-1.0, 1.0);
rsq = rsq + vec[i] * vec[i];
}
} while (rsq > 1.0);
for (i = 0; i < ndim; i++)
}

/*
* PICKBOX: pick point within box.
*/

void pickbox(real vec[], int ndim, real size)
{
int i;

for (i = 0; i < ndim; i++)
vec[i] = xrandom(- size, size);
}
```